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Development of an Innovative Algorithm 
for Aerodynamics-Structure Interaction 
Using Lattice Boltzmann Method 

EXECUTIVE SUMMARY 

The lattice Boltzmann equation (LBE) is a kinetic formulation which offers an alternative 
computational method capable of solving fluid dynamics for various systems. Major advantages 
of the method are owing to the fact that the solution for the particle distribution functions is 
explicit, easy to implement, and the algorithm is natural to parallelize. In this final report, we 
summarize the works accomplished in the past three years. Since most works have been 
published, the technical details can be found in the literature. Brief summary will be provided in 
this report. 

In this project, a second-order accurate treatment of boundary condition in the LBE 
method is developed for a curved boundary and tested successfully in various 2-D and 3-D 
configurations. To evaluate the aerodynamic force on a body in the context of LBE method, 
several force evaluation schemes have been investigated. A simple momentum exchange 
method is shown to give reliable and accurate values for the force on a body in both 2-D and 3-D 
cases. Various 3-D LBE models have been assessed in terms of efficiency, accuracy, and 
robustness. In general, accurate 3-D results can be obtained using LBE methods. The 3-D 19-bit 
model is found to be the best one among the 15-bit, 19-bit, and 27-bit LBE models. To achieve 
desired grid resolution and to accommodate the far field boundary conditions in aerodynamics 
computations, a multi-block LBE method is developed by dividing the flow field into various 
blocks each having constant lattice spacing. Substantial contribution to the LBE method is also 
made through the development of a new, generalized lattice Boltzmann equation constructed in 
the moment space in order to improve the computational stability, detailed theoretical analysis on 
the stability, dispersion, and dissipation characteristics of the LBE method, and computational 
studies of high Reynolds number flows with singular gradients. Finally, a finite difference-based 
lattice Boltzmann method is developed for inviscid compressible flows. 
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Summary 

The lattice Boltzmann equation (LBE) is an alternative kinetic method capable of solving 
hydrodynamics for various systems. Major advantages of the method are owing to the fact that 
the solution for the particle distribution functions is explicit, easy to implement, and natural to 
parallelize. Because the method often uses uniform regular Cartesian lattices in space, curved 
boundaries are often approximated by a series of stairs that leads to reduction in computational 
accuracy. In this work, a second-order accurate treatment of boundary condition in the LBE 
method is developed for a curved boundary. The proposed treatment of the curved boundaries is 
an improvement of a scheme due to Filippova & Hanel [J. Comp. Phys. 143 , 426 (1998)]. The 
proposed treatment for curved boundaries is tested against several flow problems: 2-D channel 
flows with constant and oscillating pressure gradients for which analytic solutions are known, 
flow due to an impulsively started wall, lid-driven square cavity flow, and uniform flow over a 
column of circular cylinders. The second-order accuracy is observed with solid boundary 
arbitrarily placed between lattice nodes. The proposed boundary condition has well behaved 
stability characteristics when the relaxation time is close to 'A, the zero limit of viscosity. The 
improvement can make a substantial contribution toward simulating practical fluid flow 
problems using the lattice Boltzmann method. 


Paper published in: J. Computational Physics, 155, 307-330, 1999. 
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Lattice Boltzmann Method for 3-D Flows 
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Summary 

In this work, we investigate two issues that are important to computational efficiency and 
reliability in fluid dynamic applications of the lattice Boltzmann equation (LBE): (1) 
Computational stability and accuracy of different lattice Boltzmann models and (2) the treatment 
of the boundary conditions on curved solid boundaries and their 3-D implementations. Three 
athermal 3-D LBE models (Q15D3, Q19D3, and Q27D3) are studied and compared in terms of 
efficiency, accuracy, and robustness. The boundary treatment recently developed by Filippova 
and Hanel (1998, J. Comp. Phys. 147 , 219) and Mei et al. (1999, J. Comp. Phys. 155 , 307) in 2- 
D is extended to and implemented for 3-D. The convergence, stability, and computational 
efficiency of the 3-D LBE models with the boundary treatment for curved boundaries were tested 
in simulations of four 3-D flows: (1) Fully developed flows in a square duct, (2) flow in a 3-D 
lid-driven cavity, (3) fully developed flows in a circular pipe, and (4) a uniform flow over a 
sphere. We found that while the 15-velocity 3-D (Q15D3) model is more prone to numerical 
instability and the Q27D3 is more computationally intensive, the Q19D3 model provides a 
balance between computational reliability and efficiency. Through numerical simulations, we 
demonstrated that the boundary treatment for 3-D arbitrary curved geometry has second-order 
accuracy and possesses satisfactory stability characteristics. 


Paper published in: J. Computational Phys., vol. 161, 680-699, 2000. 
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Summary 

Accurate evaluation of the hydrodynamic force on a curved body is an important issue in 
the method of lattice Boltzmann equation for fluid flow problems. This issue has not been 
systematically studied so far. The present work investigates two approaches for force evaluation: 
the method of momentum exchange and the method of stress integration. The boundary condition 
for the particle distribution functions on curved geometries is handled with second order 
accuracy based on our recent works [Mei et al. (1999, J. Comp. Phys. 155 , 307), & Mei et al. 
(2000, J. Comp. Phys. 161 , 680)]. The stress integration method is tedious in 2-D flow and 
difficult to implement in 3-D flow in general; in comparison, the momentum exchange method is 
reliable, accurate, and easy to implement in both 2-D and 3-D flows. Several test cases are 
selected to evaluate the present methods, including: 1) 2-D pressure-driven channel flow; 2) 2-D 
uniform flow over a columns of circular cylinder; 3) a channel flow over an asymmetrically 
placed circular cylinder with vortex shedding; 4) pressure-driven flow in a circular pipe; and 5) 
3-D flow over a sphere. The drag evaluated by using the momentum exchange method in the 
LBE agrees well with the exact or other published results. 


Paper currently under review for publication at Physical Review E. 
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A Multi-block Lattice Boltzmann Method for Fluid Flows 


Dazhi Yu, Renwei Mei, & Wei Shyy 
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Summary 

Compared to the Navier-Stokes equation-based approach, the method of lattice Boltzmann 
Equation (LBE) offers an alternative treatment for fluid dynamics. The LBE method often 
employs uniform lattices to maintain a compact and efficient computational procedure, which 
makes it less efficient to perform flow simulations when there is a need for high resolution near 
the body and/or there is a far-field boundary. To resolve these difficulties, a multi-block method 
is developed. In this method, the flow field is divided into blocks. In each block, the grid is 
uniform with desired resolution. In this paper, an accurate interface treatment between 
neighboring blocks is derived to ensure the continuity of mass, momentum, and stresses across 
the interface. Several test cases are employed to demonstrate that the present multi-block method 
can substantially improve the accuracy and computational efficiency. 

Paper published as AlAA 2000-2614, Fluids 2000 , Denver, Colorado, 6/19-22, 2000. 

Paper is also currently under review for publication in Int. J. of Num. Method for Fluid Flows. 
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Summary 

The generalized hydrodynamics (the wave vector dependence of the transport coefficients) of a 
generalized lattice Boltzmann equation (LBE) is studied in detail. The generalized lattice 
Boltzmann equation is constructed in moment space rather than in discrete velocity space. The 
generalized hydrodynamics of the model is obtained by solving the dispersion equation of the 
linearized LBE either analytically by using perturbation technique or numerically. The proposed 
LBE model has a maximum number of adjustable parameters for the given set of discrete 
velocities. Generalized hydrodynamics characterizes dispersion, dissipation (hyperviscosities), 
anisotropy, and lack of Galilean invariance of the model, and can be applied to select the values 
of the adjustable parameters that optimize the properties of the model. The proposed generalized 
hydrodynamic analysis also provides some insights into stability and proper initial conditions for 
LBE simulations. The stability properties of some two-dimensional LBE models are analyzed 
and compared with each other in the parameter space of the mean streaming velocity and the 
viscous relaxation time. The procedure described in this work can be applied to analyze other 
LBE models. As examples, LBE models with various interpolation schemes are analyzed. 
Numerical results on shear ow with an initially discontinuous velocity proje (shock) with or 
without a constant streaming velocity are shown to demonstrate the dispersion e_ects in the LBE 
model; the results compare favorably with our theoretical analysis. We also show that whereas 
linear analysis of the LBE evolution operator is equivalent to Chapman-Enskog analysis in the 
long- wavelength limit (wave vector k = 0), it can also provide results for large values of k . Such 
results are important for the stability and other hydrodynamic properties of the LBE method and 
cannot be obtained through Chapman-Enskog analysis. 

Paper published in Physical Review E. Vo. 61, No. 6, 6546-6562, June 2000. 
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Finite Difference-Based Lattice Boltzmann Method for Inviscid 

Compressible Flows 


Weiping Shi 1 , Wei Shyy 2 and Renwei Mei 2 

'Department of Mathematics, Jilin University, Changchun 130023, China 130023 

department of Aerospace Engineering, Mechanics & Engineering Science, University of 

Florida, Gainesville, Florida 3261 1-6250 

Summary 

A finite difference-based lattice Boltzmann model, employing the 2-D, 9-speed square (D2Q9) 
lattice for the compressible Euler equations, is presented. The model is constructed by allowing 
the particles to possess both kinetic and thermal energies. Such a lattice structure can represent 
both incompressible and compressible flow regimes. In the numerical treatment, to attain 
desirable accuracy, the total-variation-diminishing (TVD) scheme is adopted with either the 
minmod function or a second-order corrector as the flux limiter. The model can treat 
shock/expansion waves as well as contact discontinuity. Both one- and two-dimensional test 
cases are computed, and the results are compared with the exact as well as other reported 
numerical solutions, demonstrating that there is consistency between macroscopic and kinetic 
computations for the compressible flow. 

Paper is accepted for publication in Int. J. of Heat and Mass Transfer 
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The lattice Boltzmann equation (LBE) is an alternative kinetic method capable 
of solving hydrodynamics for various systems. Major advantages of the method are 
due to the fact that the solution for the particle distribution functions is explicit, easy 
to implement, and natural to parallelize. Because the method often uses uniform 
regular Cartesian lattices in space, curved boundaries are often approximated by a 
series of stairs that leads to reduction in computational accuracy. In this work, a 
second-order accurate treatment of the boundary condition in the LBE method is 
developed for a curved boundary. The proposed treatment of the curved boundaries 
is an improvement of a scheme due to 0. Filippova and D. Hanel (1998, J. Comput. 
Phys. 147, 219). The proposed treatment for curved boundaries is tested against 
several flow problems: 2-D channel flows with constant and oscillating pressure 
gradients for which analytic solutions are known, flow due to an impulsively started 
wall, lid-driven square cavity flow, and uniform flow over a column of circular 
cylinders. The second-order accuracy is observed with a solid boundary' arbitrarily 
placed between lattice nodes. The proposed boundary condition has well-behaved 
stability characteristics when the relaxation time is close to 1/2, the zero limit of 
viscosity. The improvement can make a substantial contribution toward simulating 
practical fluid flow problems using the lattice Boltzmann method. © 1999 Academic Press 


I. INTRODUCTION 


There has been a rapid progress in developing and employing the method of the lattice 
Boltzmann equation (LBE) [1-3] as an alternative computational technique for solving 
complex fluid dynamic problems (see the comprehensive reviews in [4, 5]). In a traditional 
method for computational fluid dynamics (CFD), the macroscopic variables, such as ve- 
locity u and pressure p, are obtained by solving the Navier-Stokes (NS) equations [6-8]. 
The lattice Boltzmann equation approximates the kinetic equation for the particle mass 
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distribution function f(x , £, r) on the mesoscopic level, such as the Boltzmann equation 
with the single relaxation time approximation [9], 

M + V / = -i(/-/ <0) ), (l) 

where £ is the particle velocity, / (0) is the equilibrium distribution function (the Maxwell- 
Boltzmann distribution function), and X is the relaxation time. The right hand side (RHS) 
term models the effect of the fluid viscosity on the molecular level through the collision 
(relaxation) process. The macroscopic quantities (such as mass density p and momentum 
density pu) are the hydrodynamic moments of the distribution function /, 


P = j f{x,Z,t)d 3 Z 

(2a) 

pu = J £f(x,£,t)d 3 £. 

(2b) 


It has been shown that the velocity space £ can be discretized into a finite set of points 
{£<J without affecting the conservation laws [10-12]. In the discretized velocity space the 
Boltzmann equation (1) becomes 

~ + S a -Vf a = - l -{f a -ti eq) ) (or =0,1,2 8 for 2-D) (3) 

ot A 

for the distribution function of discrete Velocities /«(*, t) = f(x , £ a , /). The equilibrium 
distribution function, f^ eq \ and the discrete velocity set {£ a } can be derived explicitly 
[ 10 — 12 ]. 

For the 2-D square lattice shown in Fig. 1, we use e a to denote the discrete velocity set, 
and we have [13] 


4 3 2 



FIG. 1. A 2-D, 9-bit (or 9-velocity) lattice. 
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e a = 0, for a = 0, 

= (cos((a — 1)7t/4), sin((a - 1)7t/4)) for a = 1, 3, 5, 7, (4) 

= V2(cos((a — 1)jt/4), sin((a — 1)jt/ 4)) for a = 2, 4, 6, 8, 


where c = 8x/8t, 8x , and <5/ are the lattice constant and the time step size, respectively, and 
= pu>„ 1 + -|e a ' w + 2^4 (*« ' “) 2 “ ^“ 2 “ “ - ( 5 ) 

where 


[ 4/9, 

a = 0 


w a = | 1/9, 

a = 1,3, 5, 7 

(6) 

i 1/36, 

a =2,4, 6,8. 


With the discretized velocity space, the hydrodynamic moments are given by 


8 

8 


11 

M 

= £/r 

(7a) 

a=0 

o=0 


and 



8 

8 


pu = E*<»/« 
0 = 1 

= E^/r. 
0=1 

(7b) 


The speed of sound of this model is c s — c/V 3, and the equation of state is that of an ideal 
gas, 

P = pc). (8) 

Equation (3) is one of numerous ways to model the transport equation of /, Eq. (1). 

Based on the Chapman-Enskog analysis, the solution for f a (x, t) may be expanded in 
the form of 

/.(*, 0 = fi eq) (x, 0 + /„ (1) (x, /) + •••, (9) 

where /j 1 * in Eq. (9) is formally smaller than in the expansion. Substitution of Eq. (9) 
into Eq. (3) leads to 

/ a (,) (jc,r) = ~x[^+e a -V/H. (10) 

ot 

Proceeding with the Chapman-Enskog analysis, it can be shown that the Euler equations 
can be recovered from the solution for and the NS equations are recovered in the 
near incompressible limit (i.e., the Mach number M = \u\/c s « 1) by the first two terms in 
Eq. (9). The viscosity of the fluid is 


v = Xc). 


( 11 ) 
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Equation (3) can be further discretized in space and time. The completely discretized 
form of Eq. (1), with the time step St and space stepe„5r, is 

fa(Xi + e a 8t , t + St) - fa(x,,t ) = ~ [fcc(Xi,t) ~ fa‘ q) (*/ , 0] , (12) 

where r = k/St, and jc ,* is a point in the discretized physical space. The above equation is 
the lattice Boltzmann equation [1-3] with Bhatnagar— Gross-Krook (BGK) approximation 
[9]. The left-hand side (LHS) of Eq. (12) is physically a streaming process for particles 
while the RHS models the collisions through relaxation. 

Although the lattice Boltzmann equation historically originates from the lattice gas cel- 
lular automata [14, 15], it is indeed a special finite difference form of the continuous 
Boltzmann kinetic equation, i.e., the LHS of Eq. (3) is discretized along the direction of the 
characteristic line with discretization of phase space and time tied together [10, 11]. The 
leading order truncation error of such a discretization is then taken into account exactly by 
modifying the viscosity in the NS equation derived from Eq. (12) to 



The positivity of the viscosity thus requires that r > 1/2. The lattice Boltzmann scheme 
consists of two computational steps, 

collision step /*(*/, t) - /*(*,, t) = — - [/ a (x,, t) — t)] (14a) 

streaming step f a (Xi +e a &t, t + St) = , f), (14b) 

where f a and f a denote the pre- and post-collision state of the distribution function, re- 
spectively. The advantages of solving the lattice Boltzmann equation over the NS equations 
can now be seen. In the kinetic equation for f a given by Eq. (3), the advection operator 
is linear in the phase space whereas the convection term is nonlinear in the NS equation. 
In traditional CFD methods, the pressure is typically obtained by solving the Poisson or 
Poisson-like equation derived from the incompressible NS equations that can be time con- 
suming. In the LBE method, the pressure is obtained through an extremely simple equation 
of state p = pc f. This is an appealing feature of the LBE method. The discretized Eq. (14) 
for f a is explicit in form, easy to implement, and natural to parallelize. The collision step 
is completely local. The streaming step takes very little computational effort at every time 
step. 

However, unlike solving the NS equations for which the non-slip condition for u on a 
solid wall is satisfied at the macroscopic level, there is no corresponding, physically based 
boundary condition for f a on a solid wall at the mesoscopic level. For a lattice node located 
on the fluid side at X/, as illustrated in Fig. 2, Eq. (14b) clearly indicates a need for the 
information of f a ztx b on the solid side. Therefore all the effort in the previous treatment 
of the boundary conditions in the LBE models is mainly focused on the calculation of f a 
moving from the wall into the fluid region. In previous works of the LBE, the most often 
used boundary condition on the wall is the so-called bounce-back scheme [16—18]. In the 
bounce-back scheme, after a particle distribution f a streams from a fluid node at x f to a 
boundary node at along the direction of e ai the particle distribution f a scatters back to 
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FIG. 2. Layout of the regularly spaced lattices and curved wall boundary. 


the node Xj along the direction of e«(=— e a ) as f a . Since the wall position x U) was forced 
to be located at jc*, this is referred to as bounce-back on the node (BBN) [19]. However, 
a finite slip velocity at the stationary wall exists [20, 18] and the accuracy for the flow 
field is thus degraded due to the inaccuracy of the boundary conditions [17]. In simulating 
suspension flows using the LBE, Ladd placed the solid walls in the middle between the 
lattice nodes [21]. This is referred to as bounce-back on the link (BBL). It has been shown 
that the BBL scheme gives a second-order accurate result for straight walls [24, 18]. Noble 
et al. developed a second-order accurate boundary condition to compute f a but it is only 
applicable to straight walls in triangular lattice space [22]. He et al. generalized the scheme 
of Noble et al to arbitrary lattice [18]. Chen et ai placed the wall on the lattice node so 
that x b is one lattice inside the wall [23]. They used an extrapolation of f a on the fluid side 
(including the wall node) to obtain f a atJC^.Zou and He proposed to apply the BBL scheme 
only for the non-equilibrium part of f a at the wall [24]. 

For a curved geometry, the use of BBL requires approximation of the curved solid bound- 
ary by a series of stair steps. The geometric integrity cannot be preserved by such an ap- 
proximation. For high Reynolds number flows, the integrity of geometry is important since 
the vorticity generation and stress distributions are sensitive to the geometrical resolution. 
To this end. He and Luo proposed to use the LBE with nonuniform grid with second order 
interpolations [10, 25, 26], He and Doolen further applied the interpolation to the LBE 
with curvilinear coordinates or body-fitted coordinates [27]. Mei and Shyy solved Eq. (3) 
in curvilinear coordinates using the finite difference method [28]. While the wall geom- 
etry is accurately preserved in body-fitted coordinates, the flexibility to handle complex 
geometries is maintained by using the numerical grid generation techniques common to 
the Navier-Stokes solvers. It should be noted that perhaps the most profound and rigor- 
ous theoretical treatment of the boundary condition along the wall is given by Ginzbourg 
and d’Humieres [29]. The scheme proposed by Ginzbourg and d’Humieres is local and 
accurate up to second order in Chapman-Enskog expansion. However, this work has not 
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attracted sufficient attention because its implementation is not as easy as the bounce-back 
scheme. 

In this work, a robust, second-order accurate treatment for the distribution function f a 
near a curved boundary is developed based on the method recently proposed by Filippova 
and Hanel (hereinafter referred to as FH) [30]. In Ref. [30], the boundary condition for 
f a on the solid side is evaluated using Eq. (3) for / a , and the Taylor series expansion in 
both space and time for f a near the wall. FH reported numerical results for a uniform flow 
over a cylinder [30]. However, it is found in this work that when tested in a pressure driven 
channel flow (see implementation and discussions in Section II) there is a strong boundary- 
condition- induced instability when the distance from the wall to the first lattice on the fluid 
side is less than half of the lattice size. 

Using the Taylor series expansion for the velocity u near the wall, a new treatment 
for f a near a curved wall is proposed in this work. While maintaining a second-order 
accuracy of the solution in handling curved walls, the computational stability is improved 
so that lower viscosity, or higher Reynolds number, can be attained in the LBE simulations. 
The new boundary condition treatment is tested systematically to assess the temporal and 
spatial accuracy and robustness in 2-D channel flow with constant and oscillating pressure 
gradients, flow due to an impulsively started wall, lid-driven square cavity flow, and flow over 
a column of circular cylinders. Detailed comparisons for the flow field are made with either 
analytic solutions or well-resolved numerical solutions of the Navier-Stokes equations by 
using a finite difference method. The improved boundary treatment represents a significant 
step towards solving practically relevant flow problems using the LBE method. 


IL FORMULATION FOR THE IMPROVED BOUNDARY CONDITION 

Filippova and Hanel [30] considered a curved boundary lying between the lattice nodes 
of spacing Sx t as illustrated in Fig. 2, and briefly presented the derivation of their scheme 
for the treatment of a curved boundary. However, they did not offer explanation to justify 
the theoretical basis of their method. It is instructive to first reexamine their derivation 
thoroughly. Based on the insight gained, an improved boundary treatment is then proposed. 


A. Reexamination of and Comments on the Filippova-Hanel Treatment 

The macroscopic flow has a characteristic length of L. The lattice nodes on the solid and 
fluid side are denoted as x b and x j, respectively, in Fig. 2. The filled small circles on the 
boundary, x Wy denote the intersections of the wall with various lattice links. The boundary 
velocity at x w > the intersection with the wall on the link between x b and Xj y is u w . The 
fraction of an intersected link in the fluid region is A, that is, 


I x f -x w \ 

I Xf -X b \ ‘ 


(15) 


Obviously, 0 < A < 1 and the horizontal or vertical distance between x b and x w is A ■ 8x 
on the square lattice. Suppose the particle momentum moving from Xj to x b is e a and the 
reversed one from x b to xj is e& = — e a . After the collision step, f a on the fluid side is 
known, but not on the solid side. (Hereafter we shall use and fa to denote the velocity 
and the distribution function coming from a solid node to a fluid node, and fa is the unknown 
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to be computed.) To finish the streaming step, 

fa( x f = x b + f + St) = f a( x b> 0* 0^) 

it is clear that fa( x b,t) is needed. To construct fa( x b> 0 based upon some known in- 
formation in the surrounding, Filippova and Hanel essentially proposed using the linear 
interpolation [30], 

fa(x b ,t) = (1 - X)fa(X / ,t) + xfa* ) (X b ,t) + 2w a pj i e & -U w , (17) 

where u w ~=u(x Wi t) is the velocity at the wall and x is the weighting factor (to be de- 
termined) that controls the linear interpolation (or extrapolation) between f a ( x f>0 an d 
f* (xb , t) y a fictitious equilibrium distribution function given by 

T 3 9 3 

0 = W a p(x f , t) 1 + - 3 e a ■ u bf + ^(e„ • “/ ) 2 - 2 ci Uf ' Uf ' 

In the above, Uj = u{xj , t) is the fluid velocity near the wall and u b f is to be chosen. It is 
emphasized here that the weighting factor x depends on how u b j is chosen. However, the 
choice of u b f is not unique. For example, either u b f =uj or a linear extrapolation using 
uyj = (A — 1)h//A + u w / A appears reasonable. 

To determine x in Eq. (17), FH considered flows under the condition 

L/(cT) « 1, (19) 

i.e., the flow has an intrinsic characteristic time scale T that is much larger than the advection 
time on the lattice scale, L/c. This “slow-flow” condition enabled FH to approximate 
fa(x f , t + St) in Eq. (16) by /«(*/, f), 

fa ( x f =x b +e a 8tyt +8t) = + + •••. 

For the purpose of the order-of-magnitude estimate, it is seen that 0(dfa/dt) = 0(fa/T) 
so that 

fa(Xf,t + &t) = fa(Xf,t) l + — fat* fit) ! + *(££) */«(*/. 0. (20) 

It is noted that under condition (19) the neglected terms are of 0(~^f) which are much 
smaller than the O(^) terms of present interest (in deriving an accurate boundary condition 
for f& (x b , t)). Applying the Chapman-Enskog expansion in the form given by Eqs. (9)— < 10) 
and invoking the “slow flow” approximation, one obtains 

f&(x f , t) = fa‘ q) (x f , t) — X O^L+ea-Vf? +■■■ 

* tf q \x f ,t) -Xe a (21) 

For given by Eq. (5), the leading order term in V fl q is given by pw a Q/c 2 )ea • Vw 
since the rest are higher order terms in the near incompressible flow limit. Noticing that 
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X = r<5f, Eq. (21) becomes 

3 

/a(x/,0 ~ fa{xj,t) - rStpw a -^e a Vu e a 

3 3 

= fa q) ( x f> 0 - 2pw a -^u f • e a - rStpu; a —e a Vu e a (22) 

which approximates the LHS of Eq. (16). To expand the RHS of Eq. (16) in terms of the 
small computational parameter 


Sx cSt 

T = ~L 


« 1, 


it is first noted that t) in Eq. (18) can be expressed as 

/j* ) ( *f>. 0 = + W a p^e a - ( u bf - Uf) 

so that the RHS of Eq. (16), or Eq. (17), can be rewritten as 


(23) 


(24) 


f a (Xb , t ) 

3 

% 0 + (1 - X)(l - 1 /T)/ a (1> Cr/,0 + w a p-^e a ■ (x«fc/ - XUf - 2 u w ) 

= fa q) ( x /' f ) - (1 “ X)( r - ^)Stw a p^e a ■ Vu e a 

3 

+ w a p— e a (x“bf ~XUf ~2u w ). (25) 

C“ 

Based on linear interpolation, Ubf ^ (A — 1)« //A -4- u w j A, expanding the velocity Uj near 
the wall (x U) ) using Taylor series, and noticing that x^—xj ~8te a , one obtains Ubf — 
Uj % Equating Eqs. (22) and (25) and matching terms linear in 8t results in 

X = (2A — l)/r. For «/,/ = u /, we have — Uf = 0 in Eq. (25). Matching to 0(8t) then 
requires x — (2 A — l)/r = (2A — l)/(r — 1). FH found that u^j — (A — l)u //A -f u w f A 
gives computationally stable results only for A > 1/2. Hence, they proposed that 

U\>f = (A — \)uj/A -fii^/A and x = (2A - 1 )/t for A > - (26) 

and 

Ubf —u f and x = (2A — 1 )/(t — 1) for A < -. (27) 

To recapitulate, there are three independent assumptions that have been made in the 
foregoing derivation. These are: (i) the Chapman-Enskog expansion in the form given by 
Eqs. (9)-(10) is valid; (ii) the intrinsic time of the unsteady flow must be large compared 
with the advection time on the lattice scale given by Eq. (19); (iii) the lattice space must 
be small compared with the characteristic length scale of the flow as given by Eq. (23) so 
that the Taylor series expansion for the velocity field near the wall is valid. There have been 
a large number of papers in the existing literature regarding the validity and usefulness of 
Chapman-Enskog expansion for the solution to the Boltzmann equation. The “slow flow” 
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k y 



FIG. 3. Lattice distribution in channel flow simulations w ith arbitrary A. 


condition is introduced to simplify the derivation of the boundary condition for f a \ the 
implication of this assumption will be briefly addressed later in comparing the computational 
results with that based on the conventional bounce-back scheme. The last assumption is a 
typical computational resolution requirement. 

Equation (17) is essentially a linear interpolation (or extrapolation) and is used continu- 
ously in the computation. When the weighting factor / becomes too large, instability may 
develop. For A > 1/2, |xl = |2A - 1 |/r is always less than 2 since the positivity of the 
viscosity in the LBE scheme requires r > 1/2. For A < 1/2, |/| = |(2A — l)/(r — 1)| and 
it may become too large when r is near 1. To illustrate this point, a fully developed pressure 
driven 2-D channel flow is considered. The grid arrangement is shown in Fig. 3. For steady 
flow, a constant pressure gradient Vp along the .x-direction is applied and can be treated as 
a body force. This is included [32] after the collision step by 

3 d 

f a (Xi , t) = fa ( Xj ,i)-w a -T -f-e a ■ x , (28) 

c l ax 

where x is the unit vector along the Jt-axis. The boundary condition for f a (x t , t) on the wall 
follows those given by Eqs. (17), (18), (26), and (27). At the inlet (i = 1) and exit (i = N x% 
in which N x is the number of lattices in the x -direction) the following is imposed, 

/«(* = 1J) = /*(*’== 2, y), (29) 

fa(i = N x J) = Mi = N x -lJ). (30) 

With Eq. (29), the velocity profile at the inlet, u x (i = 2, y), is not needed. Instead, the fully 
developed velocity profile is sought as part of the solutions. In this part of the investigation, 
M y = 35 is used. The exact solution for the velocity profile (given by Eq. (36)) is used 
as the velocity initial condition which differs from the final steady state solution due to 
numerical errors. The equilibrium distribution function f^ eq) based on the exact solution 
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FIG. 4. (a) Regions of stability and instability in the LBE computation for fully developed 2-D channel flow 

using FH’s boundary condition, Eqs. (17), (18), (26), (27), for A < 1/2. (b) Regions of stability and instability 
in the LBE computation for 2-D channel flow entrance problem using FH’s boundary condition, Eqs. (17), (18), 
(26), (27), for A <1/2. 

for the velocity profile is used as the initial condition for f a . The pressure gradient is set to 
^ = — 1.0 x 10 -6 . All computations are carried out using double precision. 

For A < 1 /2, it is found that the computation is unstable for certain range of values of 
r. Figure 4a shows the stability-instability boundaries in the (r, A) space obtained from a 
large number of computations. For A < 0.2, the computation becomes unstable when r < 1 . 
The large instability region is an apparent source of concern for FH’s scheme when A < 1 /2 
since lower viscosity can only be achieved when r is close to 1 /2. 

One may speculate that the instability in the above example results from the lack of 
specifying an inlet velocity profile, u x (y), or due to the extrapolation of f a at the inlet given 
by Eq. (29). To examine this possibility, a channel flow entrance problem is considered. 
Uniform velocity profiles, u x (y) = ~(H 2 /\2pv)(dp/dx) and //^(y) =0 in which H is the 
channel height, are specified at i = 1.5 (half-way between the first and second lattices) and 
the distribution functions f a (i = 1, j) for a = 1, 2, and 8 are obtained using Eq. (17) with 
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X =0 in accordance with A = 1/2 at i = 1.5. The boundary conditions on the wall are 
based on Eqs. ( 17), ( 1 8), (26), and (27). The exit boundary condition for the f a ’s is given by 
Eq. (30). Hence the extrapolation for f a at the inlet is completely eliminated and the velocity 
profiles at the inlet are exactly given. Two types of initial conditions are used. Whenever 
possible, the equilibrium distribution functions corresponding to the uniform inlet velocity 
are specified at t =0 throughout the flow field. This works for relatively larger values of r. 
However, instability can be encountered when r is considerably larger than the upper dash 
curve shown in Fig. 4a for the same value of A ( < 1 /2) . A second type of initial condition is 
thus implemented. A converged solution at a relatively large value of r is used as the initial 
condition for a smaller value of r. The value of r is incrementally decreased to obtain the 
converged solutions for the new, smaller values of r. When the actual instability region is 
approached, the increment in r is maintained as small as 0.01 or 0.005. In the computation, 
fa = — 1 -0 x 10~ 6 , N y = 35, and N x — 65 are used. When the Reynolds number is low (due 
to the use of the small pressure gradient and larger r), the exit velocity profile is very close 
to the exact solution corresponding to the fully developed channel flow which validates the 
solution procedure. 

The stability-instability boundaries obtained through a large number of computations 
are shown in Fig. 4b. It is noted that the stability-instability boundaries are very similar to 
that shown in Fig. 4a for the fully developed channel flow despite the dramatic difference 
in the inlet boundary condition. Thus the source of the instability must result from the 
implementation of the solid wall condition. An alternative scheme must be developed to 
overcome this shortcoming. 

B. Improved Treatment for a Curved Boundary 

We realize that the flexibility in the construction of /j #) ( 'x bt t) is the key to achieving an 
improved computational stability as well as accuracy. Since / — (2 A — l)/(r — 1) given 
by Eq. (27) leads to a larger value of x when r is close to 1, it is desirable to reduce the 
magnitude of x hy increasing the magnitude of the denominator in the expression for / . 
For A > 1/2, Ubf is the fictitious fluid velocity inside the solid and the denominator for x 
is r. For A < 1/2, u b f was chosen by FH to be uj which is the fluid velocity at x f and it 
leads to (r — 1) in the denominator for x* Thus, we propose to use Eq. (26) for A > 1/2 
and use 


u bf = Uff=Uf(x f + e & 8t,t) 


for A < 


I 

2 * 


(31) 


Thus 


Ubf ~Uf =Uf(x f + e s 8t y t)-u f (x f ,t) = StVu-e a . (32) 


This requires 


-r(l - x)(l - 1/t) X = 2A z (33) 

to match the O(St) terms in equating Eqs. (22) and (25). Hence 


X — (2A — 1)/(t — 2) for A < ~. 


(34) 
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FIG. 5. (a) Regions of stability and instability in the LBE computation for fully developed 2-D channel 
flow using the present boundary condition, Eqs. (17), (18), (26), (31), (34), for A < 1/2. (b) Regions of stability 
and instability in the LBE computation for the 2-D channel flow entrance problem using the present boundary 
condition, Eqs. (17), (18), (26), (31), (34), for A < 1/2. 


To test the improvement in the stability, the steady state, fully developed, pressure driven 
2-D channel flow is again considered. Equations (31), (34) are used in lieu of Eq. (27). The 
rest of the implementation is exactly the same as described in the last section. Figure 5a 
shows the stability-instability boundary in the (r. A) space for the fully developed channel 
flow. By comparing Fig. 5a with Fig. 4a, the improvement in the stability of the present 
treatment for this simple geometry case is clearly seen. 

For the channel flow entrance problem, boundary conditions at the inlet and exit and 
the procedure for specifying the initial conditions are the same as described in the last 
section. Equations (31), (34) are used to replace Eq. (27) for the solid wall. The stability- 
instability boundary in the (t, A) space for the entrance flow problem is shown in Fig. 5b. 
Close agreement in the stability-instability boundaries between Figs. 5a and 5b suggests 
that the improvement in the computational stability is not related to the treatment of the 
inlet boundary conditions. The improvement results rather from the different treatment in 
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the solid wall boundary condition. A direct consequence of this improvement is that lower 
values of r, or lower viscosity u, can now be used. 

One may speculate at this point that Ufix / -f- 2 t) can also be used for u^j when 
A < 1/2. This would further improve the stability since x = (2A — l)/(r — 3). This is 
correct in principle. However, since the use of u j{x f + t) as u^j already allows the 
use of i whose value is close to 1 /2, there is little practical need to use Uj that is too far 
away from the wall. 

For transient flows, a second-order extrapolation can be used for 


Ubf 


~~i ^~u(x f , t) + ^ u(x wt t) + * A [u(x w , 0 - 0 + A)||(X/, 0 

A A A(l + A) 

+ Au(xj + ez&t, /)] for A > ~ . (35) 


This treatment helps to improve the accuracy in the velocity approximation when u(x,t) is 
not well resolved near the wall. Finally, it is easily seen that the present boundary condition 
treatment can be extended to 3-D flow problems involving curved geometry. The efficacy 
of such an extension will be examined in another paper. 


III. RESULTS AND DISCUSSIONS 

For the proposed boundary condition treatment to be useful, several issues need to be 
addressed: spatial and temporal accuracy, ability to handle geometric singularity, and the 
flexibility to handle complex geometry. Channel flows with constant and sinusoidally oscil- 
lating pressure gradients with analytic solutions are used to assess the spatial and temporal 
accuracy. The Stokes first problem (i.e., the flow due to an impulsively started wall) allows 
one to examine the response of the computed flow field to an imposed singular acceleration. 
The standard lid-driven cavity flow has a bounded domain but possesses stress or vorticity 
singularities near the comers between the moving and stationary walls. Finally, flow over a 
column of circular cylinders is the case used to assess the impact of the boundary treatment 
on the accuracy of the flow field around a curved boundary. 

A. Pressure Driven Channel Flows 

At steady flow, the exact solution for the x-velocity profile is given by 

Uexact(y) = - nh (36) 

2 ax pv 

where H = N y — 3 + 2 A and r] = y/H = (j - 2 + A)///. To assess the computational error 
of the LBE solution of the velocity, u^(y ), the following relative L 2 -norm error is defined 

{Jo l u LBE(y) ”* u exact(y )] 2 dy} 

2 “ u?^u,r ' 

With the oscillating pressure gradient, ^ = Be ~ l0)t , the exact solution can be easily 
expressed in complex variables. An important parameter in this flow is the Stokes number 
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St defined as 


St = 



( 38 ) 


The Stokes number is the ratio of the channel height H to the thickness of the Stokes layer 
Jv/oj. Since the error can vary with time, a time average over one period (T =2tt/oj) is 
needed and the relative error is 

_ {Jo Jo l u LBE(y< t) ~ Ugaaiy, t)]~ dy dt} 

[/o r /o W “Lw(y.Orfy*] 1/2 

In the BGK model. At = Ax = Ay = 1. Comparing with the channel height H = N y — 3 + 
2 A, the dimensionless grid size (or grid resolution) is H l . 

Figure 6a shows the dependence of the relative L 2 -norm error on the channel height 
H for r = 0.55 and A =0.0, 0.25, and 0.5. A maximum value of N y = 131 is used. The 
second-order accuracy is demonstrated in the range of H investigated. It has been well 
established that the accuracy of the LBE method for the interior points is of second order. 
The fact that the overall accuracy is of second order in the present case means that the 
accuracy in the boundary condition is at least of second order. It is worth noticing that the 
derivation given in Section II ensures that f a is second-order accurate near the wall. It does 
not guarantee the second-order accuracy of the velocity field near the wall. To address this 
issue, the wall slip velocity, u w = u x (y = 0), is evaluated using a second-order extrapolation 
based on u x (y = A), u*(;y= 1 4* A), and u x (y= 2+ A). Since the true wall velocity in 
the pressure driven channel flow is zero, the wall slip velocity u w provides a measure of 
the accuracy for the treatment of the wall velocity. Figure 6b shows the dependence of u w , 
normalized by the centerline velocity u max = — (H 2 /Spv)(dp/dx) y on H for A =0.0, 0.25, 
and 0.5 with r = 0.55. Quadratic convergence is clearly observed in all three cases which 
demonstrates the second-order accuracy of the velocity field near the solid wall. This is 
entirely consistent with the results shown in Fig. 6a which involves global convergence 
rather than the local (y = 0) convergence. Figure 6c shows the relative error as a function 
of A using the present boundary treatment (Eqs. (17), (18), (26), (31), (34)) for 0 < A < 1. 
The error in the range of 0 < A < 1/2 is comparable to that in the range of 1/2 < A < 1. 
The present boundary condition treatment does not induce larger computational error and 
is substantially more robust. Furthermore the second-order accuracy is achieved in general 
by the present treatment for A < 1/2. 

Figure 7 shows the dependence of the relative L 2 -norm error on the channel height H 
in the oscillating pressure driven channel flow for Stokes number Sr = Hy/cojv = 1 and 8. 
For St — 1, the Stokes layer is as thick as the channel height H. For A = 0.25, 0.5, and 0.75, 
second-order accuracy in space is clearly demonstrated. Since the time step 8t in LBE is 
equal to the spatial resolution 8x , the accuracy in time must also be of second-order in 
order for the time-averaged L 2 -norm error to have a slope of 2 in Fig. 7. For St — 8, the 
Stokes layer thickness is about 1/8 of the channel height so that the computational error 
due to the insufficient resolution of the Stokes layer is a significant part of the error. For 
A = 0.25, the first lattice in the flow field is only a quarter of the lattice size away from 
the wall. The Stokes layer is thus better resolved for A = 0.25 (denoted by solid circles in 
Fig. 7) than for A = 0.5 and 0.75. However, as H increases, the difference between A = 0.25 
and A = 0.5 and 0.75 becomes smaller since all have reasonable resolutions in the Stokes 





FIG. 6. (a) Dependence of relative L 2 -norm error on the lattice resolution H = N y — 3 + 2A, in steady state 
pressure-driven channel flow simulations, (b) Quadratic convergence of the wall slip velocity in steady state 
pressure-driven channel flow simulations, (c) Relative L 2 -noTm error as a function of A in steady state pressure- 
driven channel flow simulations. 


321 




322 


MEl, LUO, AND SHYY 



Channel height H in lattice unit 


FIG. 7. Dependence of the L 2 -norm error on the lattice resolution H — N y — 3 -F 2A in oscillating pressure 
driven channel flow. Stokes number St — H +J7oJv. 


layer. Although the slope for the error curve for A = 0.25 is observed to be about 1.5 that 
is less than 2, it is an indication of the better- than-expected accuracy at the low resolu- 
tion end. 

B. Stokes First Problem: Flow Due to an Impulsively Started Wall 

Fora wall located at y = 0 that is impulsively started, an unsteady Stokes layer of thickness 
0(*/vi) develops near the wall. For a fixed-grid computation, the error at small time is 
expected to be large due to insufficient spatial resolution. In the LBE method, this is also 
compounded by the use of fixed 8t (=<$* = <$y = 1). Figure 8 shows the velocity profiles 
at t = 100 (in lattice unit). The wall velocity is V =0.1 in lattice unit. The relaxation time 
r =0.52 gives kinematic viscosity v = 0.0067. Similar to the oscillating pressure driven 
channel flow, the error is smaller for A = 0.25 than for A =0.5 and 0.75 due to a better 
spatial resolution near the wall. Figure 9 shows the temporal variation of the relative L 2 - 
norm error defined as 


_ t)] 2 dy} l/2 0 

for A =0.25, 0.5, and 0.75. The result using the standard bounce-back on the link (BBL) 
scheme, which always sets A = 0.5, is also shown. The large relative errors in the beginning 
are due to the smaller values of the denominator in the above equation. It should be empha- 
sized that this flow at small time is difficult to deal with for any computational technique due 
to the singular acceleration and large spatial gradient. For an impulsively started Couette 
flow, the long-time solution approaches the exact linear velocity profile because the LBE 
method is a second-order accurate one. It is interesting to note that the present boundary 
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FIG. 8. Velocity profiles at r = 1 00 (lattice unit) of an impulsively started plate with various values of A. The 
bounce back on the link (BBL) always sets A = 1/2. 


condition treatment for A = 0.5 gives a slightly smaller error than the BBL scheme in this 
highly transient case. In such a transient flow, the computational accuracy in the near-wall 
region is typically dictated by the near-wall spatial resolution which must be smaller than 
the Stokes layer thickness in order to resolve the local flow field. In a finite difference 
calculation for such a flow, St and Sx can be independently chosen. If 8x is not sufficiently 
small, further reduction in St will not lead to improvement in accuracy. At small f, neither 
the BBL scheme nor the present treatment resolved the Stokes layer so that the error is large. 
After the Stokes layer grows to a certain thickness, the spatial resolution becomes adequate 



FIG. 9. Relative L 2 -norm error of the velocity profile u x (y ) during the initial transient of the impulsively 
started plate with various values of A. The “linear” version of the boundary condition corresponds to Eq. (26). 
The “quadratic” version corresponds to Eq. (35). The BBL is limited to A =0.5 only. 
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FIG. 10. Velocity profiles at the center ( x/H — 1/2) in lid-driven cavity flow with various values of A at 
Re = 100. 


and the accuracy then improves. In view of the “slow flow” condition (19) introduced in 
the derivation, the performance of the current boundary treatment is comparable or better 
than the conventional bounce-back on the link scheme. 


C. Flow in a Lid-Driven Square Cavity 

Figure 10 shows the velocity profiles at the center (x/H = 1/2) of the cavity of width H 
at Re = 100 with r = 0.6. Only 35 x 35 lattices are used and the cavity width is H — N x — 3-b 
2 A = 32 -f- A. This requires the lid velocity to be V = vRe/H = 3.33 /H in the lattice unit. It 
has a negligible compressibility effect for H ~ 32. A well-resolved finite difference solution 
for the velocity field based on the stream function- vorticity formulation is also shown for 
comparison. The velocity profile with A = 0. 1 agrees well with the finite difference solution. 
For A = 0.5, the result is rather reasonable with such a resolution. The difference is slightly 
larger on the negative velocity part for A = 0.9. The corner singularity in stress (or vorticity) 
is well handled for r — 0.6 and N x = 35. However, for t close to 0.5 and with N x = 35, 
the corner singularity induces wiggles in the velocity field. This issue will be examined 
elsewhere. The flow field for Re= 1000 is obtained with 67 x 67 lattices using A = 0.1, 
0.5, and 0.9. Similar behavior in the velocity profiles is observed. 


D. Uniform Flow over a Column of Circular Cylinders 

To simulate the external flow over a single cylinder would require placing the outer 
boundary far away from the cylinder. In order to keep the computational effort at a reason- 
able level in using constant space lattices, a column of circular cylinders of radius r and 
center-to-center distance H is considered instead. The flow field that needs to be computed 
is thus limited to —H < y < H. At y = — //, the lattice is j = 2. The boundary conditions 
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at j = i for the M s are given by the following symmetry properties, 

/o0\ 1) = Mi, 3), /,(/, 1) = Mi, 3), Mi, D = Mi, 3), 

Mi, 1) = Mi, 3), Mi, 1) = Mi, 3), /s(u D - Mi, 3), (41) 

/fiO\ 1) = Mi, 3), / 7 (|\ 1) - Mi, 3), /g(i\ 1) = Mi, 3). 

Similar conditions hold at y — H for j — N y . At the inlet, the uniform velocity, u = V, is 
specified at i — 1.5. Using A = 0.5, / = 0, Eq. (17) is applied to obtain the f a 9 s at i = 1. At 
the exit, a simple extrapolation is used, 


f a (N Xy j) = 2 f a (N x - 1, j) - f a (N x - 2, j) for a - 4, 5, and 6. (42) 


On the surface of the circular cylinder, Eqs. (17), (18), (26), (31), and (34) proposed in this 
paper are used to update the boundary conditions for the Ms. 

Figure 11 shows the velocity profile u(x = 0, y)/V for H/r =20at/?e = 2Vr/v = 10 
using r = 3.5. Two values of relaxation time r (=0.505 and 0.525) are used. For r = 3.5, there 
are only 7 lattices from the front to the back stagnation points. The finite difference solution 
is obtained using body-fitted coordinates [33] and over 200 grid points are distributed along 
the upper surface of the circle. These two solution with r = 0.505 and r = 0.525 are virtually 
identical to each other and they are both close to the finite difference solution. Figure 12 
shows the centerline (y = 0) velocity variations, upstream and downstream, respectively, 
at Re= 10 and 40. The sharp gradient near the front stagnation point, the length of the 
separation bubble, the maximum of the separation bubble velocity, and the recovery of 
the wake velocity are all in excellent agreement with the well resolved finite difference 
solution. 

As can be seen now, an important improvement of the present boundary condition treat- 
ment over the bounce-back scheme is that it can preserve the accuracy of the geometry 



FIG. 11 . Velocity profiles at x ~ 0 for uniform flow over a column of cylinders. The cylinder has a diameter 
(2r) of 7 lattice units. The cylinder center-to-center distance H = 70 lattice units. 
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FIG. 12. Centerline ( y =0) velocity variation for a uniform flow over a column of cylinders. Finite difference 
results are based on co — if formulation and are well resolved, (a) Upstream; (b) downstream. 


under consideration. To further demonstrate this point, consider flow over a circular cylin- 
der of radius (r) with the coordinate centered at the center of the cylinder. For r = 3.4 
and 3.8, the front stagnation points are located at x = —3.4 and 3.8, respectively. With 
the bounce-back on the link (BBL) scheme, the front stagnation points in both cases will 
be placed atx = -3.5 which is half-way between the lattice at x = — 4 and x = -3 on the 
centerline. In the present method, A = 0.6 and 0.2 for r = 3.4 and 3.8, respectively. The 
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FIG. 13. Comparison of the velocity profiles at x = 0 for r = 3.0, 3.2, 3.4, 3.5, 3.6, 3.8, and 4.0 for Re — 1 0 
and H /r = 20. 


difference in A can be accurately incorporated in the evaluation of fa(x b , t). This implies 
that although the boundary links for r = 3.4 will be different from those for r — 3.8, the 
flow fields based on r — 3.4 and r = 3.8 should be nearly the same when the coordinates 
are normalized by the radius r. To validate this point, a series of computations are carried 
out for r = 3.0, 3.2, 3.4, 3.5, 3.6, 3.8, and 4.0 for H/r =20 at fle = 10. The profiles of 
the dimensionless x-component velocity u x /U as a function of y/r at x = 0 are compared 
for these seven different radii r in Fig. 13. Excellent agreement is observed. Figure 14 
compares the u x /U as a function of x/r at y = 0 for both the downstream and upstream 
regions for these seven different radii. Again, all seven cases compare very well even in the 
near wall region. This clearly demonstrates that the present boundary condition treatment 
has maintained geometric fidelity even with coarse grid resolutions. 

It is noted that the interpolation for fa(x b , t) given by Eq. (17) is carried along the line 
in the direction of e a . The results for flow over a cylinder are quite satisfactory. Other 
interpolation procedures can certainly be devised to use more information on neighboring 
lattices in the flow field. However, this will result in a lot more complications in the imple- 
mentation. It is not clear if such an attempt will necessarily lead to further improvement 
over the present approach. 


IV. CONCLUSION 

In this work a second-order accurate boundary condition treatment for the lattice 
Boltzmann equation is proposed. A series of studies are conducted to systematically vali- 
date the accuracy and examine the robustness of the proposed boundary condition in steady 
and unsteady flows involving flat and curved walls. Compared with the existing method for 
treating boundary condition in the lattice Boltzmann method, the proposed treatment has 
the following advantages: (i) It can preserve the geometry of interest without truncating 
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a 



b 



FIG. 14. Comparison of the centerline velocity for r = 3.0. 3.2, 3.4, 3.5, 3.6, 3.8, and 4.0 for Re — 10 and 
H/r = 20. (a) Upstream region; (b) downstream region. 


it into a series of stair steps, (ii) The boundary treatment generally results in solutions of 
second-order accuracy for the velocity field in space, and in time for some cases, (iii) Com- 
pared with the widely used bounce-back on the link scheme, the present treatment gives 
comparable or better results for the flow field under otherwise identical computational 
parameters. 
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In this work, we investigate two issues that are important to computational effi- 
ciency and reliability in fluid dynamic applications of the lattice Boltzmann equation 
(LBE): (1) Computational stability and accuracy of different lattice Boltzmann mod- 
els and (2) the. treatment of the boundary conditions on curved solid boundaries and 
their 3-D implementations. Three athermal 3-D LBE models (Q15D3, Q19D3, and 
Q27D3) are studied and compared in terms of efficiency, accuracy, and robustness. 
The boundary treatment recently developed by Filippova and Hanel (1998, J. Comp . 
Phys. 147, 219) and Mei el al (1999, J. Comp . Phys. 155, 307) in 2-D is extended to 
and implemented for 3-D. The convergence, stability, and computational efficiency 
of the 3-D LBE models with the boundary treatment for curved boundaries were 
tested in simulations of four 3-D flows: (1) Fully developed flows in a square duct, 
(2) flow in a 3-D lid-driven cavity, (3) fully developed flows in a circular pipe, and 
(4) a uniform flow over a sphere. We found that while the 15- velocity 3-D (Q15D3) 
model is more prone to numerical instability and the Q27D3 is more computationally 
intensive, the Q19D3 model provides a balance between computational reliability 
and efficiency. Through numerical simulations, we demonstrated that the boundary 
treatment for 3-D arbitrary curved geometry has second-order accuracy and possesses 
satisfactory stability characteristics. © 2000 Academic press 

Key Words: lattice Boltzmann equation; boundary condition for curved geometries; 
accuracy; 3-D flows. 


I. INTRODUCTION 

1.1. Basic Notion of the Lattice Boltzmann Equation 

In one fashion or another, conventional methods of computat ional fluid dynamics (CFD) 
compute pertinent flow fields, such as velocity u and pressure p, by numerically solving the 
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Navier-Stokes equations in space* and time t [1-3]. In contrast, various kinetic methods 
use the transport equation, or the Boltzmann equation in particular, for various problems in 
fluid dynamics. The Boltzmann equation deals with the single particle distribution function 
/(*, £, /), where £ is the particle velocity, in phase space (x y £) and time t. Recently, the 
method of the lattice Boltzmann equation (LBE) [4-7] has become an alternative to the 
conventional CFD methods employing Navier-Stokes equations. The theoretical premises 
of the LBE method are that (1) hydrodynamics is insensitive to the details of microscopic 
physics, and (2) hydrodynamics can be preserved so long as the conservation laws and 
associated symmetries are respected in the microscopic or mesoscopic level. Therefore, 
the computational advantages of the LBE method are attained by drastically reducing the 
particle velocity space £ to only a very few discrete points without seriously degrading 
hydrodynamics. This is possible because the LBE method rigorously preserves the hydro- 
dynamic moments of the distribution function /, such as mass density and momentum 
fluxes, and the necessary symmetries [8-10]. 

One popular kinetic model is the Boltzmann equation with the single relaxation time 
approximation [11], 


M + £. V / = -i[/-/<°>], (l) 

where f is the particle velocity, / (0) is the equilibrium distribution function (the Maxwell- 
Boltzmann distribution function), and k is the relaxation time. The mass density p and 
momentum density pu are the first (D - f 1) hydrodynamic moments of the distribution 
function / and / <0 \ where D is the dimension of velocity space. 

To solve for / numerically, Eq. (1) is first discretized in the velocity space f using a finite 
set of velocities {£*} without affecting the conserved hydrodynamic moments [9-11], 

^ + la ' V/ a = - i [/„ - /„«“«>] . (2) 


In the above equation, / a (jc, /) = /(*, t) and /j eq) = / (0) (*, £ a , t) are the distribution 
function and the equilibrium distribution function of the otth discrete velocity respec- 
tively. The 9- velocity (or 9-bit) LBE model on the 2-D square lattice, denoted as the Q9D2 
model, has been widely used for simulating 2-D flows. For 3-D flows, there are several 
cubic lattice models, such as the 15-bit (Q15D3), 19-bit (Q19D3), and 27-bit (Q27D3) 
models [12], which have been used in the literature. All three models have a rest particle 
(with zero velocity) in the discretized velocity set {£«}. A minor variation of those models 
is to remove the rest particles from the discrete velocity set; the resulting models are known 
as the Q14D3, Q18D3, and Q26D3 models, respectively. The LBE models with a rest 
particle generally have better computational stability. For athermal fluids, the equilibrium 
distributions for the Q9D2, Q15D3, Q19D3, and Q27D3 models are all of the form [9] 


/i eq) = pw a 


3 9 3 

1 + —e a u + — (e„ • u ) 2 - —u 


2c 4 


2c 2 


(3) 


where w a is a weighting factor and e a is a discrete velocity, c = Sx /St is the lattice speed, and 
Sx and St are the lattice constant and the time step, respectively. (The values of the weighting 
factor w a for the Q15D3, Q19D3, and Q27D3 models and the diagrams illustrating the lattice 
structures for the Q 15D3 and Q 19D3 models are given in the Appendix.) It can be shown that 
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fW is in fact a Taylor series expansion of the Maxwellian / <0) [8, 9]. This approximation 
of / (0) by the above f^ cq) makes the method valid only in the incompressible limit u/c -*■ 0. 

With the velocity space discretized, the hydrodynamic moments of / and / (0) are eval- 
uated by the following quadrature formulas: 

= (4a) 
a a 

pu = ^e 0 / 0 = E e “^o <eq) - 

a a 


The speed of sound of the above 3-D LBE models is c s = c/%/3 and the equation of state is 
that of an ideal gas p — pc f. The viscosity of the fluid is v = kc*. 

Equation (2) is often discretized in space, a:, and time, r, into 

f a (Xi +e a 8t,t +8t) - fa(XiJ) = -i[/a(*n 0 " (*/ , 0] . ( 5 ) 

where r -k/St. This is the lattice Boltzmann equation with the Bhatnagar-Gross-Krook 
(BGK) approximation [1 1] and is often referred to as the LBGK model [4, 5]. The viscosity 
in the NS equation derived from Eq. (5) is 

v = (r — l/2)c^<$f. (6) 

This choice for the viscosity makes the LBGK scheme formally a second order method 
for solving incompressible flows [9]. The positivity of the viscosity requires that r > 1/2. 
Equation (6) can be solved in the following two steps: 

collision step: / a (*n 0 = /(*(*/, t) — “[/<*(*! , 0 — /j eq) » 0] > (^ a ) 

streaming step: f a (Xi +e a 8t, t + <5f) = f a ( x ii 0. (?b) 

where f a denotes the post-collision state of the distribution function. It is noted that the 
collision step is completely local, and the streaming step is uniform and requires little 
computational effort. Equation (7) is explicit, easy to implement, and straightforward to 
parallelize. 

1.2. Boundary Condition on a Solid Surface 

To date, most Neumann-type boundary conditions for a solid boundary used in the LBE 
method are based upon the bounce-back boundary condition: A particle colliding with a sta- 
tionary wall simply reverses its momentum. Much of the previous work on LBE boundary 
conditions is devoted to the analysis and improvement of the bounce-back boundary condi- 
tion [13-21, 27]. The bounce-back boundary condition can attain second-order accuracy if 
the boundary is fictictiously placed halfway between two nodes. That is, the second-order 
accuracy of the bounce-back boundary condition can only be achieved when the boundaries 
are located right in the middle of two neighboring lattices [A = 0.5; see Eq. (8)]. (Readers 
are referred to our recent work [22] for a summary of the previous work.) This prevents 
the direct application of the bounce-back-type boundary conditions to simulate a solid 
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FIG. 1. 2-D projection of the layout of the regularly spaced lattices and curved wall boundary. 


body with smooth curvature. To circumvent this difficulty, Mei and Shyy solved Eq. (2) 
in curvilinear coordinates using a finite difference method to solve for f a [28]. One can 
also use body-fitted curvilinear coordinates with interpolation throughout the entire mesh, 
except at the boundaries where the bounce-back boundary condition is used [29]. In more 
recent works [22, 23], Cartesian coordinates are adopted with interpolation used only at the 
boundaries. These techniques rely on the freedom of using interpolation techniques. We 
used the latter technique in the present work. 

As shown in Fig. 1 for a 2-D projection involving a 3-D body, the streaming step requires 
the knowledge of /«(**,, f). in which e s = —e a , at x b on the solid side in order to obtain 
fa(x/, St) for the lattice node located on the fluid side atjc / = x b + e^St . Defining 


I x f -x w \ 

\Xf-X b \ 


( 8 ) 


as the fraction of an intersected link in the fluid region, it is seen that 0 < A < 1 and the 
horizontal or vertical distance between x b and x w is ( 1 — A )Sx on the cubic lattice. 

Based on the work of Filippova and Hanel [23], hereinafter referred to as FH, Mei et ai 
[22] proposed the following treatment for f s ( x by t) on curved boundaries: 

3 

fa(x b ,t) = (1 - x)f a (Xf,t) + xfa\ x b<t) + 2w a p—e a ■ u w (9) 


with 


f«\xb,t) = U) a p(Xf,t) 


3 9 3 

1 + -? e a ■ u bf + ^(e a - Uf ) 2 - ^-fUf ■ Uf 


( 10 ) 


and 


u b f = (A - l)u fj A + Mu,/ A and x = (2A - 1 )/z 

u bf = Ufj = Uf(x f + e s St,t) and x = (2 A - l)/(r - 2) 


for A >1/2 (11) 

for A < 1 /2. (12) 
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It is noted that Eq. (12) for u bf and x differs from that originally proposed by FH. The 
choice for u b f given by Eq.( 12) improves the computational stability for r < 1 and A < 1/2 
[22], Since Eqs. (9M 12) are in vector form, they can be directly extended to 3-D flows with 
curved boundaries. 

1 .3. Scope of the Present Work 

The present study examines two issues in 3-D incompressible fluid dynamics simulations 
with arbitrary boundaries using the LBE method: (i) The performance of various 3-D 
athermal LBE models for viscous flows, and (ii) the efficacy and reliability of the extension of 
the curved boundary treatment from 2-D to 3-D flows. We focus on the stability and accuracy 
of the computation and the robustness in handling an arbitrary curved geometry. In Section II, 
a modification of the choice of u bf and the expression for / when A > 1/2 is proposed in 
order to further improve the computational stability of the boundary treatment. In Section III, 
numerical results for four 3-D steady flows are examined and various computational issues 
are addressed. These four cases are: (i) pressure driven fully developed flow in a square 
duct; (ii) 3-D lid-driven cavity flow; (iii) pressure driven fully developed flow in a circular 
pipe; and (iv) uniform flow over a sphere. In cases (i) and (iii), the LBE-based numerical 
solutions can be compared with known exact solutions so that the accuracy of the LBE 
solutions can be determined. The difference in these two cases is that A is a constant in the 
square duct while A varies around the solid boundary in the circular pipe. In the lid-driven 
cavity flow, the singularities at corners between the moving and stationary walls allows for 
a performance assessment of various LBE schemes. The flow past a sphere is an external 
flow around a 3-D blunt body. In all four cases, detailed assessments are made in terms of 
error norms and velocity profiles. It will be demonstrated that accurate and robust solutions 
are obtained using the newly proposed boundary conditions along with the selected LBE 
models. 


II. MODIFICATION OF THE BOUNDARY CONDITION FOR A > 1/2 

Equations (9)-( 12) are first applied to a fully developed pressure driven 2-D channel flow 
by using the 3-D LBE model Q19D3. At the inlet (i = 1) and exit (/ = N x , in which N x is 
the number of lattices in the x-direction) the following zero derivative condition is imposed 
after the collision step: 


fad = 1 ,j,k) - fa d = 2 ,j,k), (13) 

/„( i = N X J, k) = f a (i = N X - 1,7. k). (14) 

At k = 1 and k = N z , the same is imposed: 

fad, j, k = l) = fad, j, k = 2), (15) 

fad, j, k = N z ) = fad, j, k = N z — 1). (16) 


The constant pressure gradient V p along the x-direction is treated as a body force and is 
included in the solution procedure after the collision step and the enforcement of the above 
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FIG. 2. Stability boundary of FHs scheme in a square duct flow for A near 1 . 


zero-derivative conditions as 

3 dp 

fa (X| » 0 == fa (X| i 0 'X, (17) 

c~ ax 

where x is the unit vector along the x-axis. On the solid walls (v =0 and y = H ), Eqs. (9)- 
(12) are used. The exact solution for the velocity is used as the velocity initial condition. The 
equilibrium distribution function /j eq) based on the exact solution for the velocity profile 
is used as the initial condition for f a . The pressure gradient is set to ^ — 1.0 x 10 -6 . All 

computations are carried out using double precision. 

It was found that the computations are stable for r close to 0.5 (for example, r = 0.505) 
as long as A is not too close to unity (for example, A < 0.87). When A is equal to 1, stable 
computation can only be carried out for r no smaller than 0.6. Figure 2 shows the stability- 
instability boundary for the channel flow simulation with a system size N x x N y x N z =5x 
35 x 5, near A = 1. Also shown by the dashed line is the stability-instability boundary for 
the channel flow simulation using the Q9D2 model and with a system size N x x N y = 5 x 35, 
near A = 1. It is clear that similar behavior exists in both 2-D and 3-D channel flow simu- 
lations. When the computation for the pressure driven flow in a square duct was carried out 
using the Q19D3 formulation, a similar stability-instability boundary was encountered. 

Ideally, one would like to use a fixed value of r for the entire range of 0 < A < 1 in a 
simulation. Computational stability would then require the use of r around 0.6, instead of 
a value that is close to 0.5, which makes it difficult to simulate a lower viscosity or higher 
Reynolds number flow. To overcome the restriction imposed by the numerical stability 
requirement due to interpolation, it would be useful if one could decrease the value of 
X = (2A — l)/r given by Eq. (11). This can be accomplished by using 

u bf = [1 — 3/(2 A)]a/ + 3/(2A)u w and X = (2 A - l)/(r + 1/2) for A > 1/2. 

(18) 

That is, the velocity u^f is evaluated at (x/, 4* 1/2 e a ), instead of atx/,, using the information 
atxy andXw through linear extrapolation. 
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With Eq. (18) replacing Eq. (1 1), the channel flow simulations using the Q19D3 lattice 
model are carried out again for A from 0.85 to 1. Satisfactory results for the velocity pro- 
files are obtained for r = 0.505 with N x x N y x N z = 5 x 35 x 5 in terms of computational 
stability. For A < 0.85, the accuracy of the solutions using Eqs. (11) and (18) is the same 
when the computations are stable. 


III. RESULTS AND DISCUSSIONS 


3.1. Fully Developed Flow in a Square Duct 

For fully developed flow inside a square duct of height// defined by the region -a<y<a 
and — a < z ^ u, where a — H /2, the axial velocity profile can be found in Ref. [24, p. 123]. 


16 a 2 dp ^ I" cosh(rt7Tz/2a) 

u x (y>z)- 3 dx 2-, cosh(/nr/2) 

n=1.3 t 5,... 

Figure 3 compares the exact axial velocity profiles at z = 0 and the LBE-based solution 
with A = 0.2 and H = 2a = 32.4. A total of N x x N y x N z = 13 x 35 x 35 grid points are 
used. The pressure gradient is = —1.0 x 10 -6 and r =0.52. The 19-bit model is used 
in the simulations. Excellent agreement was obtained. 

Figure 4a shows the dependence of relative L 2 -norm error, 


cos(tt7ry/2fl) 


(19) 


E 2 = 


{lo Jo [ w lbeCv> z) - w exact (y i z)Y dy dz^ 

r ZZ Z _ . 1 l / 2 


1/2 


( 20 ) 


fo"fo H “exact (y,z)dydz^ 
on the duct height or the lattice resolution H = N y - 3 + 2A. The integral is evaluated by 



FIG. 3. Comparison of axial velocity profiles in a pressure driven square duct flow at z = 0 between the exact 
solution and the LBE-based solution with A = 0.2. 
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FIG. 4. (a) Dependence of relative L 2 -norm error on the lattice resolution H — N y - 3 + 2A in steady state 

pressure-driven duct flow simulations, (b) Relative Ln-norm error as a function of A in steady-state pressure-driven 
duct and channel flow simulations. 

the trapezoidal rule. As was demonstrated by Mei etal [22], the boundary treatment results 
in second order convergence for 2-D channel flow. Figure 4a clearly shows that the total 
error (from both the flow field and the boundary condition) of the LBE solution in 3-D flow 
decays quadratically. 

Figure 4b shows the relative L 2 -norm error E 2 as a function of A in the duct flow using 
13 x 35 x 35 grid points and r = 0.52. For the purpose of comparison, the relative L 2 -norm 
error in the 2-D channel flow simulation using the Q9D2 model with N y = 35 and r = 0.52 
is also shown. The relative error is larger in 3-D duct flow than in the 2-D channel flow. 
Nevertheless, the error exhibits the same qualitative behavior in both 2-D and 3-D as a 
function of A. 
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It should be noted that the accuracy of the Q9D2 model and the Q19D3 model is different 
in the sense that beyond the conserved moments (density and momentum in athermal fluids), 
these two models have different accuracy in preserving higher order moments (fluxes) 
[9, 10]. The Q9D2 model preserves all the moments up to second order in momentum 
space, which include momentum fluxes, and maintains the isotropy of these moments, 
whereas the QI9D3 model can preserve density and momentum, but cannot maintain the 
same accuracy and isotropy of the fluxes like the Q9D2 model does. The only 3-D equivalent 
of the Q9D2 model in terms of accuracy of the moments is the Q27D3 model [9, 10]. 


3.2. Simulation Results for 3-D Lid-Driven Cavity Flows 

Lid-driven cavity flow has been studied extensively in the CFD community. Most research 
has been focused on 2-D problems. Limited numbers of reliable numerical results for 
steady state 3-D cavity flows have been obtained in the past several years. In this study, the 
multiblock finite difference solution of the NS equations obtained recently by Salom [25] 
is used to compare with the present LBE-based results. 

The size of the cavity is W \ the number of grids is N x x N y x A-, and N x = N y = N z . 
The driving lid is placed at y = H, moving along the direction of x-axis with a speed 
U =0.1 in lattice units. Figure 5a compares profiles of horizontal velocity u x (y) obtained 
using 33 x 33 x 33 lattices with the solution to the NS equations at x/H -z/H = 0.5 for 
Re = 400. All three LBE models (15-, 19-, and 27-bit) are used. For the 15-bit model, 
the computation becomes unstable and blows up at this Reynolds number with 33 3 lattice 
resolution and A = 0.5. For A = 0.5, the 19- and the 27-bit models give very similar u x (y) 
profiles and both underpredict slightly the magnitude of the minimum in the profiles. The 
19-bit model is also used with A = 0.25; there is a slight overshoot in the velocity profiles 
in comparison to the results in Ref. [25]. Figure 5b compares u x (y) profiles obtained using 
the 15- and 19-bit lattice models on the 67 3 lattice grids and A =0.5 with the NS solution 
[25] at x/H =z/H = 0.5 for Re = 400. Excellent agreement is observed. Clearly, the 19-bit 
model is superior to the 15-bit model. Although the 15-bit model requires 21% less CPU 
time and storage than the 19-bit model per lattice, it is not as robust as the 19-bit model 
and may actually require more CPU time and memory to obtain a reasonable solution since 
more lattice points are clearly needed. 

It should be noted that the stability property of the 19- and the 15-bit models is signifi- 
cantly different. All LBE models have inherent spurious invariants because of their simple 
dynamics [30]. However, the stability of the LBE models, which is very much affected 
by these spurious invariants, differs from one model to another and also depends on other 
factors such as boundary conditions and the local Reynolds number [30]. Among the three 
3-D LBE models (Q15D3.Q 1 9D3, and D27D3), the Q15D3 model is the least isotropic and 
therefore is more prone to numerical instability. This is independently verified in a recent 
work by Kandhai et al. [26]. It was observed that the Q15D3 model may induce artificial 
checkerboard invariants which are the eigenmodes of the linearized LBGK collision oper- 
ator at wave vector k = tt; this can cause spatial oscillations to develop in the flow field at 
high Reynolds number [30], Although it was pointed out that the presence of solid walls can 
suppress the oscillation in certain cases, the solid walls in the present case actually excite 
the oscillation by producing shear stress singularities at the corners between the moving and 
the stationary walls. Clearly, the Q19D3 model is better suited to handle flow singularities 
than the Q15D3 model in this case. 
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FIG. 5. Comparison of u x profiles using (a) 33 x 33 x 33 and (b) 67 x 67 x 67 lattices with a Navier-Stokes 
(NS) solution at xjH =zjH = 0.5 for to = 400 in a lid-driven cavity flow. 

Figure 6a compares the profiles of transversal velocity u y (x) obtained from various 
3-D LBE models using 33 3 lattices (grids) with the NS solution at y/H = z/H =0.5 for 
Re = 400. For A = 0.5, we found that the results from the 27-bit model deviate more from 
the NS results of Ref. [25] than the results of the 19-bit model with the same resolution in 
the spatial region 0.1 <x/H < 0.6. Both models underpredict the extrema of the velocity 
profile compared to the NS solution of Ref. [25]. For A = 0.25, the results of the 19-bit 
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FIG. 6. Comparison of u y profiles using (a) 33 x 33 x 33 and (b) 67 x 67 x 67 lattices with an NS solution 
at y/H =zj H =0.5 for tfe = 400 in a lid-driven cavity flow. 


model slightly overpredict the extrema, also shown in Fig. 6a. However, the difference is 
relatively smaller in both cases. Figure 6b shows velocity profiles with a resolution of 67 3 
grid points and the same Reynolds number Re = 400. With 67 3 lattice resolution, the result 
of the 15-bit model significantly differs from the results of the 19-bit model and that of the 
NS solution in Ref. [25]. These comparisons further suggest that the 19-bit model is better 
than the 15-bit model in terms of accuracy and stability and better than the 27-bit in terms 
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FIG. 7. Effect of Reynolds number on the centerline (xjH — z/H — 1/2) velocity profiles, u x /U , based on 
the 19-bit LBE solutions in a lid-driven cavity flow. 


of computational efficiency. The 19-bit model represents a good compromise in terms of 
both computational efficiency and reliability. 

Figures 7 and 8 show the effect of Reynolds number (from 100 to 2000) on the pro- 
files of horizontal velocity u x (y ) at x/H =z/H =0.5 and transversal velocity u y (x) at 
y j H = z/ H =0.5 based on the Q19D3 model. For Re — 100, 400, and 1000, A =0.5 is 
used. It is worth noting that for Re — 2000, the system size of 67*' , U =0.1, and r = 0.50325, 
the LBE simulation with A = 0.5 eventually becomes unstable, although the steady-state 
result of Re — 1000 is used as the initial condition for Re = 2000. When A = 0.25 is used on 
the 67 3 lattice system, no computational instability occurs and the steady-state solution is 
obtained. Weak spatial oscillations in the u x (y) andwy(x) velocity profiles were observed for 
Re = 2000, which indicates that further increase inRc would require better spatial resolution. 
It is also worth pointing out that when FH's boundary condition [23] is used for Re — 2000 
with A = 0.25, the solution eventually blows up even when converged results (based on the 
present boundary condition for A = 0.25) at Re = 2000 are used as the initial condition. 

3.3. Fully Developed Flows inside a Circular Pipe 

Figure 9 shows the 2-D projection of the discretized domain and the boundary nodes x b 
(denoted by solid symbols) on the yz plane for flow inside a circular pipe of radius R =9.5 
lattice units. Geometrically, the LBE simulation of the pipe flow differs from that of the duct 
flow in that the fraction of the intersected link A is not constant over the entire boundary. 
As seen in Fig. 4b, computational error can vary with A in the duct flow and the difference 
in the error can easily be as large as a factor of four for 0 < A < 1. Furthermore, the error is 
the smallest when A is between 0.3 to 0.6. Hence, it is reasonable to expect that the overall 
error in the solution will depend on the distribution of A in the entire set of A. 
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FIG. 8. Effect of Reynolds number on the y-component velocity profiles at yj H = z/H = 1/2 based on the 
19-bit LBE solutions in a lid-driven cavity flow. 



FIG. 9. Boundary nodes \ b (solid symbols) for flow in a pipe of radius 9.5 lattice units. 



LATTICE BOLTZMANN METHOD 


693 



FIG. 10 . Variation of relative errors in the velocity profiles as a function of pipe radius. 


Figure 10 shows the relative Z^-norm error for the axial velocity profile defined as 

_ {E(vy, J *)en[«LBE(>’;.Zi)-«exact(3' J ,Zt)] 2 } 1/2 

E 2 = £ -j J 2 , (2 1) 

^ exact 7 » 2 *)] 

where Q is the set of the discrete lattice grids inside the pipe, as a function of radius R 
for R = 3.5, 4.5, 5.5, 9.5, 13.5, 18.5, and 23.5. The pressure gradient is % = -1.0 x 10~ 6 
and r = 0.52. It is noted that each simple summation in Eq. (21) is slightly less than the 
exact integration over the entire circle due to the discretization. To ensure that such a 
treatment does not affect the qualitative behavior of the error measurement, the centerline 
axial velocity, w r , is also compared with the exact solution and the error is defined as: 


E c 


! w c,LBE u c , exact I 

I u c , exact I 


( 22 ) 


It is seen that E c behaves very similarly to E 2 and both are nonmonotonic. This oscillatory 
behavior could be due to the difference in the distribution of A, which in turn results in the 
difference of the dissipation due to the interpolation around the boundary. Shown also in 
Fig. 10 is the error E 2 of the square duct flow solution (with A = 0.2) as a function of equiv- 
alent radius H/ 7T 1/2 , which exhibits a quadratic convergence. Despite the nonmonotonic 
behavior, it can still be seen that on average, E 2 and E c decay quadratically with increasing 
radius and the accuracy in the pipe flow simulation is comparable to that in the square duct 
flow simulation. 

Figure 1 1 shows the axial velocity profiles in the pipe for R = 3.5, 5.5, 9.5, and 13.5 in 
comparison with the exact solution. Even for a very small radius R = 3.5, the LBE solution 
agrees with the exact solution remarkably well. A noticeable discrepancy in the velocity 
profile at R =9.5 is also observed in E 2 and E c shown in Fig. 10. 
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u x 

FIG, 11. Comparison of the axial velocity profiles between the LBE-based and the exact solutions for flow 
inside a cicular pipe. 

3.4. Simulation Results for a Uniform Flow over a Sphere 

The conventional LBE scheme uses uniform meshes. Without local mesh refinement, it 
is difficult to compute the external flow over a blunt body efficiently since a large number 
of grid points in the far field will be wasted. As a first attempt, the flow over a sphere is 
computed within a finite region in the transversal directions. 

As shown in Fig. 12, the outer boundary is placed at y = ±H/2 and z — ztH/ 2. At 
y~-Hf 2, the lattice is j =2. The boundary conditions at j — 1 for f a 's are given by the 
following linear extrapolation: 

fad, 1 , *) = 2 /*(!. 2, k) - fa (I, 3, *). (23) 

The velocity at j = 2 is set as 


«(/, 2, k) = w(c 3, k). 


(24) 



FIG. 12. Schematic for uniform flow over a sphere. 
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FIG. 13. Comparison of the velocity profiles at x = 0 for r = 3.0, 3.2, 3.4, 3.6, 3.8, and 4.0 for Re — 10 and 
H/r = 5. 


Similar treatment is applied at y= H /2 and z = ±H/2. The extrapolation condition given 
by Eqs. (23) and (24) allow the flow to leave the outer boundary. This helps to reduce the 
effect of the outer boundary on the flow field and on the drag force. At the inlet, a uniform 
velocity profile is imposed at ( = 1.5 (halfway between the first and second lattice points) 
and Eq. (9) is applied to obtain the condition for f a ( 1 , j , k) with x = 0. At the exit, a simple 
extrapolation is used: 


f a (N x J, k ) = 2f a (N x - 1, j, k) - f a (N x - 2, j , k). (25) 

On the surface of the sphere, Eqs. (9), (10), (12), and (18) proposed in this work are used 
to update the boundary conditions for f a ’s. Only the 19-bit LBE model is used to simulate 
the flow over a sphere. 

Figure 1 3 shows the velocity profile u x (y) based on a series of computations carried out for 
several values of the radius R =3.0, 3.2, 3.4, 3.6, 3.8, and 4.0 for H / R = 10 at Re — 10. The 
results are obtained with r = 0.7. Figure 14 compares the axial velocity profile (at y = z = 0) 
for the same set of parameters. It is worth noting that the present LBE computation does 
not have sufficient resolution for the given Reynolds number. Yet the velocity profiles agree 
with each other accurately. The fact that we have obtained a spatially accurate solution 
over a range of radii strongly suggests that the present boundary condition treatment for 
curved geometry in the LBE method is capable of handling more complex geometries while 
maintaining good accuracy. 


VI. CONCLUDING REMARKS 

Three 3-D LBE models, including the 15-, 19-, and 27-bit models, have been assessed in 
terms of efficiency, accuracy, and robustness in lid-driven cavity flow. While accurate 3-D 
results can be obtained by using various LBE models, the 19-bit model is found to be the 
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FIG. 14. Comparison of the centerline velocity (y = 0) for r = 3.0, 3.2, 3.4, 3.6, 3.8, and 4.0 for Re — 10 and 
U jr = 5. 

best for the cases investigated. The 15-bit model exhibits velocity oscillations and is prone 
to computational instability. The more complicated 27-bit model does not necessarily give 
more accurate results than the 19-bit model with the same spatial resolution. 

In this study, we have also modified the boundary condition treatment for the LBE method 
proposed by Filippova and Hanel [23] and Mei et al. [22] when the fraction of the intersected 
link on the boundary A is greater than one half. This improves the computational stability 
when A is close to 1 and r close to 1/2. 

The simulations for flows in a square duct and in a circular pipe indicate that the current 
boundary condition treatment for curved geometries results in second-order accuracy in 
3-D flows. The velocity profiles for flow over a sphere show good self-consistency of the 
solution over a range of sphere radii used. 


APPENDIX 


The Q15D3 model has the following set of discrete velocities: 


€ a 


(0, 0, 0), a = 0; rest particle 

(±1,0, 0)c, (0, ±1, 0)c, (0, 0, ±l)c, a = 1, 2, . . . , 6; group I (Al) 

(±1, ±1, ±l)c, a = 7, 8 14; group III 


and the weighting factor w a is [12] 


(2/9, 
w a = < 1/9, 

1 1/72, 


or = 0; rest particle 
a = 1, 2, . . . , 6; group I 
a = 7, 8, . . . , 14; group III. 


(A2) 


The Q19D3 model has the following set of discrete velocities: 
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FIG. Al. Discrete velocity vectors for the QI5D3 and Q19D3 lattices. 




( 0 , 0 , 0 ), 

(±1,0, 0)c, (0, ±1, 0)c, (0, 0, ±l)c, 

(±1, ±1, 0)c, (±1, 0, ±l)c, (0, ±1, ±l)c, 


a = 0, rest particle 
a = 1,2, . . . , 6; group I (A3) 
a = 7, 8, . . . , 18; group II 


and the weighting factor w a is [9] 


r 1/3, 

w a = \ 1/18, 
l 1/36, 


a = 0; rest particle 
a = 1, 2, .... 6; group I 
a = 7, 8, .... 18; group II. 


(A4) 


The Q27D3 model has the following discrete velocities: 

a = 0; rest particle 
a = 1,2, .... 6; group I 
a = 7, 8, .... 18; group II 
a = 19, 20, ... , 26; group III 

(A5) 


(( 0 , 0 , 0 ), 

I (±1,0, 0)c, (0, ±1, 0)c, (0, 0, ±l)c, 

I (± 1 , ± 1 , 0)c, (± 1 , 0, ± l)c, (0, ± 1 , ± 1 )c, 
[(±l,±l,±l)c, 
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w a — 
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, is [9] 



8/27, 

a = 0, rest particle 

2/27, 

a = 1,2,... 

, 6; group I 

1/54, 

a = 7,8,... 

, 18; group II 

1/216, 

a = 19,20, . 

. . , 26; group 


(A6) 


In the above, c = 8x /St, 8x and 8t are the lattice constant and the time step size, respectively. 
The lattice structures for the Q15D3 and Q19D3 models are shown in Fig. Al. 
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Abstract 

The present work investigates two approaches for force evaluation in the lattice Boltzmann equa- 
tion: the momentum exchange method and the stress integration method on the surface of a 
body. The boundary condition for the particle distribution functions on curved geometries is han- 
dled with second order accuracy based on our recent works [1, 2]. The stress integration method 
is computationally laborious for two-dimensional flows and in general difficult to implement for 
three-dimensional flows, while the momentum exchange method is reliable, accurate, and easy to 
implement for both two-dimensional and three-dimensional flov/s. Several test cases are selected 
to evaluate the present methods, including: (i) two-dimensional pressure-driven channel flow; (u) 
two-dimensional uniform flow past a column of cylinders; (iii) two-dimensional flow past a cylinder 
asymmetrically placed in a channel (with vortex shedding); (iv) three-dimensional pressure-driven 
flow in a circular pipe; and (v) three-dimensional flow past a sphere. The drag evaluated by using 
the momentum exchange method agrees well with the exact or other published results. 
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I. INTRODUCTION 


A Background of the lattice Boltzmann equation method 

The method of lattice Boltzmann equation (LBE) solves the microscopic kinetic equation 
for particle distribution function /(x,£,t), where £ is the particle velocity, in phase space 
(x, £) and time t , from which the macroscopic quantities (flow mass density p and velocity 
u) are obtained through moment integration of /(x,£,i). Because the solution procedure 
is explicit, easy to implement and parallelize, the LBE method has increasingly become an 
attractive alternative computational method for solving fluid dynamics problems in various 
systems [3-6]. The most widely used lattice Boltzmann equation [3-6] is a discretized 
version of the model Boltzmann equation with a single relaxation time approximation due 
to Bhatnagar, Gross, and Krook (BGK model) [7] 

3J + £-V/ = I[/-/ (0) ], (1) 

where / (0) is the Maxweil-Boltzmann equilibrium distribution function and A is the relax- 
ation time. The mass density p and momentum density pu are the first (D+ 1) hydrodynamic 
moments of the distribution function / and in D-dimensions. It can be shown that the 
particle velocity space £ can be discretized and reduced to a very small set of discrete ve- 
locities {£ Q |a; = 1,2,... ,6}, and the hydrodynamic moments of f and as well as their 
fluxes can be preserved exactly, because the moment integral can be replaced by quadrature 
exactly up to a certain order in £ [8-11]. With velocity space £ properly discretized, Eq. (1) 
reduces to a discrete velocity model (DVM) of the Boltzmann equation 

dtf a +Z a -Vf a = j{fa-fi 0) ], (2) 

In the above equation, f a (x y t) = and fa\x,t) = ont) are the distri- 

bution function and the equilibrium distribution function of the a-th discrete velocity £ Q? 
respectively. Equation (2) is then discretized in space x and time t into 

fa{ x i + “ fa( x ii 0 = 0 ” /q 0] i (3) 


2 


where r = A }5 t is the dimensionless relaxation time and e Q is a discrete velocity vector. 
The coherent discretization of space and time is done in such a way that Sx = e Q St is 
always the displacement vector from a lattice site to one of its neighboring sites. The 
equilibrium distribution function fa q \x l , t ) in the lattice Boltzmann equation (3) is obtained 
by expanding the Maxwell-Boltzmann distribution function in Taylor series of u up to second 
order [8, 9], and can be expressed in general as 



( 4 ) 


where c = 6 x /5 t ] S x is the lattice constant of the underlying lattice space; and coefficient 
w a depends on the discrete velocity set {e Q } in D spatial dimensions. In what follows, 
we shall use the lattice units of S x — 1 and S t = 1. Appendix A provides the details of 
coefficient w a and the discrete velocity set {e Q } for the two-dimensional nine- velocity model 
(D2Q9) and the three-dimensional nineteen- velocity model (D3Q19) [12], Figure 1 shows 
the discrete velocity sets of the two models. It should be pointed out that there exist other 
discrete velocity sets {e Q } which have the sufficient symmetry for the hydrodynamics [8, 9], 
A comparative study of three three-dimensional LBE models including the fifteen-velocity 
model (D3Q15), the nineteen-velocity model (D3Q19), and the twenty-seven- velocity model 
(D3Q27), in terms of accuracy and computational efficiency has been conducted by Mei et 
al. [2], It was found that the nineteen-velocity model (D3Q19) offers a better combination 
of computational stability and accuracy. The D2Q9 and D3Q19 models will be used in 
this study for force evaluation in two-dimensional and three-dimensional flows, respectively. 
Equation (3) is conveniently solved in two steps 

collision: f a {xi,t) = f a {xi,t ) - - [ f a {x u t ) - f£ q) (xi,t)] , (5a) 

r 

streaming: / a (x* -1- e a S u t + S t ) = / a (xj, t ) . (5b) 


which is known as the LBGK scheme [3, 4]. The collision step is completely local, and 
the streaming step is uniform and requires little computational effort, which makes Eq. (5) 
ideal for parallel implementation. The simplicity and compact nature of the LBGK scheme, 
however, necessitate the use of the square lattices of constant spacing (5 X = 5 y ), and con- 
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sequently lead to the unity of the local Courant-Friedrichs-Lewy (CFL) number, because 

S t = S x . 

B Boundary condition for a curved geometry in the LBE method 

Consider a part of an arbitrary curved wall geometry, as shown in Fig. 2, where the filled 
small circles on the boundary, x w , denote the intersections of the boundary with various 
lattice- to-lattice links. The fraction of an intersected link in the fluid region, A, is defined 
by 

_ 11^/ ~ X w|l (0) 

II*/ - ®&ll ' 

Obviously the horizontal or vertical distance between Xf, and x w is A S x on the square lattice, 
and 0 < A < 1. In Eq. (5b), the value of f a (x u t) needs to be constructed according to the 
location of the boundary and the boundary conditions, if the grid point Xj = Xb lies beyond 
the boundary. In the past, the bounce-back boundary condition has been use to deal with a 
solid boundary in order to approximate the no-slip boundary condition at the solid boundary 
[13-21], However, it is well understood that this bounce-back boundary condition satisfies 
the no-slip boundary condition with a second-order accuracy (for Couette and Poiseuille 
flows) at the location one half lattice spacing (A = 1/2) outside of a boundary node where 
the bounce-back collision takes place; and this is only true with simple boundaries of straight 
line parallel to the lattice grid. For a curved geometry, simply placing the boundary half-way 
between two nodes will alter the geometry on the grid level and degrade the accuracy of the 
flow field and the force on the body at finite and higher Reynolds number. To circumvent 
this difficulty, Mei and Shyy solved Eq. (2) in curvilinear coordinates using a finite difference 
method to compute f Q [22]. He and Doolen used body-fitted curvilinear coordinates with 
interpolation throughout the entire mesh, except at the boundaries where the bounce-back 
boundary condition is used [23]. In the recent works of Filippova and Hanel [24] and Mei et 
al. [1, 2], a second-order accurate boundary condition for curved geometry was developed 
in conjunction with the use of Cartesian grid in order to retain the advantages of the LBE 
method. An interpolation scheme is employed only at the boundaries to obtain f a {xi,t). 
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The detailed assessment on the impact of the boundary condition on the accuracy of the 
flow field have been given in Ref. [1] for some two-dimensional flows and in Ref. [2] for some 
three-dimensional flows. 

C Force evaluation and related works 

In spite of numerous improvement for the LBE method during the last several years, one 
important issue that has not been systematically studied is the accurate determination of 
the fluid dynamic force involving curved boundaries. Needless to say, accurate evaluation of 
the force is crucial to the study of fluid dynamics, especially in fluid-structure interaction. 
Several force evaluation schemes, including momentum exchange [15, 17] and integration of 
surface stress [23, 25], have been used to evaluate the fluid dynamic force on a curved body 
in the context of the LBE method. 

He and Doolen [23] evaluated the force by integrating the total stress on the surface of the 
cylinder and the components of the stress tensor were obtained by taking respective velocity 
gradients. Even though the body-fitted grid was used, an extrapolation was needed to obtain 
the stress in order to correct the half-grid effect due to the bounce-back boundary condition. 
Filippova and Hanel [24] developed a second-order accurate boundary condition for curved 
boundaries. However, the fluid dynamics force on a circular cylinder asymmetrically placed 
in a two-dimensional channel was obtained by integrating the pressure and deviatoric stresses 
on the surface of the cylinder by extrapolating from the nearby Cartesian grids to the 
solid boundary [24, 25]. To gain insight into the method of surface stress integration, it 
is instructive to examine the variation of the pressure on the surface of a circular cylinder 
at finite Reynolds number obtained by using the LBE method for flow over a column of 
cylinders (see Ref. [1], and Sec. Ill B) . Figure 3 shows the pressure coefficient 

^ _ P ~ Poo 

° F = W r 

on the surface obtained by using second order extrapolation where is the far upstream 
pressure. Only those boundary points, x w , intersected by the horizontal or vertical velocities, 
i.e., ei, e 3 , e 5 , and e 7 , are considered in the result given by Fig. 3. If the boundary points 
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intersected by the links in the diagonal velocities, i.e., e-z, e$, and e%, are also considered, 
the variation of C P would be more noisy. The components of the deviatoric stress tensor 
show a similar noisy pattern. It is not clear how the noise in the pressure and stresses 
affect the accuracy of the fluid dynamic force in the stress integration method. While the 
programming in the extrapolation and integration is manageable in two-dimensional cases, 
it is rather laborious in three-dimensional cases. 

Instead of the stress integration method, Ladd used the momentum exchange method to 
compute the fluid force on a sphere in suspension flow [15]. In the flow simulation using 
the bounce-back boundary condition, the body is effectively replaced by a series of stairs. 
Each segment on the surface has an area of unity for a cubic lattice. The force on each 
link [half-way between two lattices at Xf and x\, — ( Xf + e a 5t) in which x (, resides in the 
solid region] results from the momentum exchange (per unit time) between two opposing 
directions of the neighboring lattices 

T *[®a/a(®/) e & f & (Xf -\- Gadf)] 

Ot 

in which e a = — e Q . Whereas the momentum exchange method is very easy to implement 
computationally, its applicability and accuracy for a curved boundary have not been system- 
atically studied. To recapitulate, there are two major problems associated with the method 
of surface stress integration. First, the components of stress tensor are often noisy on a 
curved surface due to limited resolution near the body and the use of Cartesian grids. The 
accuracy of such a method has not been addressed in the literature. Second, the implemen- 
tation of the extrapolation for Cartesian components of the stress tensor to the boundary 
surface and the integration of the stresses on the surface of a three-dimensional geometry 
are very laborious in comparison with the intrinsic simplicity of the lattice Boltzmann simu- 
lations for flow field. The problems associated with the method of the momentum exchange 
are as follows, (a) The scheme was proposed for the case with A = 1/2 at every boundary 
intersection x w . Whether this scheme can be applied to the cases where A ^ 1/2, when, for 
example, the boundary is not straight, needs to be investigated, (b) As in the case of stress 
integration method, the resolution near a solid body is often limited and the near wall flow 
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variables can be noisy. If one uses the momentum exchange method to compute the total 
force, it is not clear what the adequate resolution is to obtain reliable fluid dynamic force 
on a bluff body at a given (moderate) value of Reynolds number, say, Re « O(10 2 ). 

D Scope of the present work 

In what follows, two methods for the force evaluation, i.e., the stress integration and the 
momentum exchange methods, will be described in detail. The shear and normal stresses on 
the wall in a pressure driven channel flow will be first examined to assess the suitability of the 
momentum exchange method when A ^ 1/2 and analyze the errors incurred. The results on 
the drag force for flow over a column of circular cylinders using these two methods will be 
subsequently assessed for the consistency. The drag coefficient at Re = 100 will be compared 
with the result of Fornberg [26] obtained by using a second-order accurate finite difference 
scheme with sufficient grid resolution. For flow over a cylinder asymmetrically placed in a 
channel at Re = 100, the unsteady drag and lift coefficients were computed and compared 
with the results in the literature. The momentum exchange method is further evaluated for 
three-dimensional fully developed pipe flow and for a uniform flow over an two-dimensional 
array of spheres at finite Reynolds number. We found that the simple momentum exchange 
method for force evaluation gives fairly reliable results for the two-dimensional and three- 
dimensional flows. 

II. METHODS FOR FORCE EVALUATION IN LBE METHOD 
A Second-order accurate no-slip boundary condition for curved geometry 

The analysis of boundary conditions for a curved boundary in the lattice Boltzmann 
equation is accomplished by applying Chapman-Enskog expansion for the distribution func- 
tion at the boundary. The following approximation for post-collision distribution function 
on the right-hand-side of Eq. (5b) can lead to a second-order accurate no-slip boundary 
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condition [1, 2, 24] 


3 

h{x b ,t) = (1 - X)Ia{Xf,t) + Xf*a(Xb,t) + 2 w a p- 2 e & ■ u w , (7) 

where 


r 3 

fa(Xb,t) = W a p(Xf,t) 1 + ^( e Q 

\ 9 . ^2 ^2 

■ U b f ) + ^( e a-Uf) 


3 

= fS^(xj, t) + W a p(xj, t)-^e a ■ ( u bf - u f ) , 

(8) 

Ubf = Uff = Uf(x f + e^St^t) , 


(9a) 

1 3 

u bf = “ Quf + ’ 


(9b) 


The above treatment is applicable for both the two-dimensional and three-dimensional lattice 
Boltzmann models. 

By substitution of Eq. (8), Eq. (7) becomes 

fa[xb,t ) = fc{x f ,t) - x [/a (*/>*) ~ fa q) ( x f^) 

3 

+w a p{xf,t)—e a -{u b f-Uf-2u w ). (10) 

Thus, the above treatment of curved boundary can be thought as a modification of the 
relaxation (the viscous effect) near the wall (via parameter x)> i n additional to a forcing 
term accounting for the momentum exchange effect due to the wall. 

B Force evaluation based on stress integration 

He and Doolen [23] evaluated the force by integrating the total stresses on the boundary 
of the cylinder dQ: 

F= f dAh-{-p\ + pv[(V:u) + {V:u) T }} , (11) 

Jd n 

where fi is the unit out normal vector of the boundary dfh In Ref. [23], a body fitted 
coordinate system together with grid stretching was used such that a large number of grids 
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can be placed near the body to yield reliable velocity gradient diUj. In general, since u is not 
the primary variable in the LBE simulations and the evaluation of u using ^ a Ca/a based 
on f a 's suffers the loss of accuracy due to the cancellation of two close numbers in f a s the 
evaluation of the derivative diUj will result in further degradation of the accuracy. Filippova 
[25] used similar integration scheme to obtain the dynamic force on the body for the force 
on a circular cylinder [24] except that the deviatoric stresses were evaluated using the non- 
equilibrium part of the particle distribution function [see Eq. (13) below]. However, since 
the Cartesian grid was used, the stress vectors on the surface of the body (with arbitrary 
A) have to be computed through an extrapolation procedure based upon the information in 
the flow field. This leads to further loss of accuracy for finite lattice size S x when the shear 
layer near the wall is not sufficiently resolved. 

In Eq. (11), the pressure p can be easily evaluated using the equation of state p = c 2 p. 
For D2Q9 and D3Q19 models, c z s = 1/3 so that p = p/3. The deviatoric stress for two- 
dimensional incompressible flow 

Tij = pv ( diUj + djUi) (12) 

can be evaluated using the non-equilibrium part of the distribution function = [f a — 

fi eq) ] 

Tij = (l - ^ ^2 /a neq) (®^) (e Q)f e QJ - ~e a ■ eJi^j , (13) 

where e a ti and e a j are ith and jth Cartesian component of the discrete velocity e a , respec- 
tively. For the flow past a circular cylinder, a separate set of surface points on the cylinder 
can be introduced in order to carry out the numerical integration given by Eq. (11). The 
values of the pressure and each of the six components of the symmetric deviatoric stress 
tensor on the surface points can be obtained using a second-order extrapolation scheme 
based on the values of p and at the neighboring fluid lattices. The force exerting on the 
boundary dfl is computed as 

F = [ dAfi- { —p\ + pu[{V:u) + (V :u) T ]} extrapolated • (14) 

Jd n 
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It is worth commenting here that for the two-dimensional flow past a cylinder, nearly half 
of the entire code was taken up by the above force evaluation procedure. 


C Method based on the momentum exchange 


In order to employ the momentum exchange method efficiently, two scalar arrays, w(i,j) 
and w b (i,j) are introduced. A value of 0 is assigned to w(i,j) for the lattice site (i,j) that 
are occupied by fluid; a value of 1 is assigned to w{i,j) for those lattice nodes inside the 
solid body. The array w b (i,j) is set to zero everywhere except for those boundary nodes, 
x b , where a value of 1 is assigned. For a given nonzero velocity e Q , e & denotes the velocity 
in opposite direction, i.e. , = — e Q (see Fig. 2). For a given boundary node x b inside the 

solid region with w b (i,j) = 1 and = 1 , the momentum exchange with all possible 

neighboring fluid nodes over a time step St = 1 is 


a^O 


fai y *Eb 0 


[1 - w(x b + e a S t )] . 


Simply summing the contribution over all boundary nodes x b belonging to the body, the 
total force (acted by the solid body on the fluid) is obtained as 


all ajfc a^O 


fa(Xb') 0 "t” fa 


[1 -iu(x6 + e„5 t )]. 


(15) 


In the momentum exchange method the force F is evaluated after the collision step is 
carried out and the value of f & at boundary given by Eq. (7) has been evaluated. The 
momentum exchange occurs during the subsequent streaming step when /g(a3f,,t + St) and 
f a (xf,t + S t ) move to x/ and x b , respectively. As mentioned in the introductory section, the 
effect of variable A is not explicitly included, but it is implicitly taken into account in the 
determination of fa{ x bi t + St). The applicability of Eq. (15) will be examined and validated. 

Clearly, the force is proportional to the number of boundary nodes x b in the above formula 
of F and the number of the boundary nodes increase linearly with the size of the body in 
a two-dimensional flow. However, since the force is normalized by pU 2 r in the formula for 
Cd in two-dimensions [see Eq. (24)], the drag coefficient Co should be independent of r. 
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III. RESULTS AND DISCUSSIONS 


For straight walls, there is no doubt that Eq. (11) together with the equation of state 
for pressure and Eq. (13) for T{j gives accurate result for the force provided that f a s are 
accurately computed. To demonstrate the correctness of Eq. (15) based on the momentum 
exchange for an arbitrary A, we first consider the pressure driven channel flow (see Fig. 4) for 
which exact solutions for the velocity and stresses are known. The second case considered is 
the two-dimensional flow past a column of circular cylinders at Reynolds number Re — 100 
and H/r = 20, where H is the distance between the centers of two adjacent cylinders. The 
values of the drag computed using the two force evaluation methods are then compared with 
the result of Fornberg [26]. The dependence of the drag on the radius r in the momentum 
exchange method is examined to assess the reliability of this method. The third case is the 
two-dimensional flow over a circular cylinder that is asymmetrically placed in a channel at 
Re = 100 (with vortex shedding). The time dependence of the drag and lift coefficients is 
compared with results in literature. 

We also consider two cases of three-dimensional flow. The first case is the pressure driven 
flow in a circular pipe for which the exact solutions for both the velocity profile and wall 
shear stresses are known. The assessment for the momentum exchange method for three- 
dimensional flows will be made first in this case. Finally, the momentum exchange method 
will be evaluated by considering the drag on a sphere due to a uniform flow over a sphere 
in a finite domain. The details for the flow field computation can be found in Ref. [1, 2]. 

A Two-dimensional pressure-driven channel flow 

In the case of the channel flow, the force on the top wall (y = H) at a given location x 
{i = jV x /2 + l, say) can be evaluated using the momentum exchange method as follows. The 
wall is located between j — N y and N y — 1 (Fig. 4). The x and y components of the force 
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on the fluid at the top wall near the ith node are 

F x = [/«(*, j) + h(i -1,3 ~ 1)] e 6,x + [fs(i,3) + h (i + 1 >i - 1)] e 8 ,x (16a) 

Fy = + Mi -1,3- 1)] e 6 ,y + [h{i,j) + h{i + 1,3 ~ 1)] e 8 , y 

+m, 3) + Mi, 3 -l)]e7,y, ( 16b ) 

where e a j denotes the jth Cartesian component of velocity e a . Since 5 X = 1, F x and F y are, 
effectively, the total shear and normal stresses, cr xy and <7 yy , which include the pressure and 
the deviatoric stresses, on the fluid element at y = H. 

Based on Eq. (13), the deviatoric component of the fluid shear stresses at j = N y — 1 
(or y = jV y - 3 + A) and N y - 2 (or y = N y - 4 + A) can be exactly evaluated based on 
the non-equilibrium part of the distribution functions in the flow field if they are correctly 
given. A linear extrapolation of the deviatoric shear stresses to y = H = N y — 3 + 2 A yields 

T = r xy (j = N y - 1) + A [r xy (j = N y - 1) - r xy (j = N y - 2)] , (17) 

where the superscript “(neq)” denotes the value computed from fa q \ the subscript w refers 
to the value at the wall. The deviatoric normal stress, Tyy!i2\ can be similarly computed. 
In a fully developed channel flow, the normal component of the deviatoric stress r yy (y) is 
expected to be zero while the total normal stress cr yy (y ) is equal to the negative pressure 
(— p). It needs to be pointed out that this method of evaluating Tx™w given by Eq. (17) for 
two-dimensional channel flow is equivalent to the method of the surface stress integration 
based on the extrapolated pressure and the deviatoric stresses on the solid wall except that 
no numerical integration on the solid surface is needed. 

After the velocity profile u x (y) is obtained from / a , the shear stress r xy on the wall can 

also be calculated using the near wall velocity profile as, 

du x _ (2 + A) [0 - u x (j = N y - 1)] 

PU dy y=H - pU ( 1 + A) A 

~P v Tf = N v - !) “ u *(3 = Ny- 2)1 • ( 18 ) 

In the above, a linear extrapolation is employed to evaluate the velocity derivative y=H 
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at the wall. Finally, the exact solution for the fluid shear stress on the wall (y = H) is 

exact = \ d V H # = I\L — 3 + 2A (19) 

'xy,w 2 dx V 

based on the parabolic velocity profile or simple control volume analysis. This exact result 
can be used to assess the accuracy of the aforementioned methods for the force evaluation. 

In the LBE simulations, the pressure gradient is enforced through the addition of an 
equivalent body force after the collision step [2, 23]. While the velocity field given by the 
LBE solution can be unique, the pressure field [thus the density field p(x, y )] can only be 
unique up to an arbitrary constant. In view of Eq. (18), it is difficult to compare the stresses 
for different cases if p{i,j) converges to different values in each case. To circumvent this 
difficulty, the density field in the channel flow simulation is normalized by p(i = 2, j = N y j 2) 
at every time step. This normalization procedure results in p(x, y) = 1 throughout the entire 
computational domain. It is also applied to the three-dimensional flow in a circular pipe. 

Table I compares the numerical values of the shear stress for a typical case (N y = 35, 
dp/dx = -10“ 6 in the lattice units, and r - 0.6) based on: given by Eq. (19), F x 

given by Eq. (16a), r xy e $ given by Eq. (17), and pv^-\ y=H given by Eq. (18). Also listed 
is the comparison between F y given by Eq. (16b) and — p. All computations are carried out 
with double precision accuracy. 

It is noted that r is identical to for a11 values of A- Closer examination of the 

shear stress profile using Eq. (13) across the channel reveals that T xy e t w (y) is also equal to the 
exact shear stress profile T®* act (y), which is linear, despite the errors in the velocity profile 
u x (y) for all values of A. A linear extrapolation, Eq. (17), for a linear profile therefore gives 
the exact wall shear stress. Thus, the exactness of r xy *w in the LBE simulation of channel 
flow indicates the reliability of the LBE solution for the stress field r/" eq) (x, y) by using 
Eq. (13). However, as Fig. 3 indicates, the accuracy of the integrating r^ neq) (x,?/) to obtain 
the fluid dynamic force in nontrivial geometries is not clear; this will be further investigated 
in the following sections. 

For 0 < A < 1, the normal force F y given by Eq. (16b) based on the momentum exchange 
method agrees exactly with the pressure on the wall. This is a rather special quantity 
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since deviatoric component of the force is identically zero. Nevertheless, the method of the 
momentum exchange does give a reliable value for the normal stress. 

For the shear (tangential) force, it is observed from Table I that for fixed dp/dx , F x 
does not change as A increases from 0.01 to 0.99. On the other hand, the exact result 
exact = 1 &(N - 3 + 2A), increases linearly with A. Further computations were carried 

xy,w 2 dx v V n 

out over a range of N y (= 35, 67, 99, and 131) and r (= 0.505, 0.51, 0.52, 0.6, 0.7, 0.8, 0.9, 
1.0, 1.2, 1.4, and 1.6). The results indicate that the momentum exchange method gives the 
shear stress on the top wall as 


1 dp 

2 dx 



(20) 


That is, F x is independent of r and A. The error in F x is zero when A = 1/3. The absolute 
error attains the maximum when A = 1 which gives the relative error of ^ for F x . Although 
the frequently used momentum exchange method is a natural choice for the force evaluation 
in conjunction with the bounce-back boundary condition for A = 1/2, one must be aware 
of that this method is not exact and the error in the force evaluation using the momentum 
exchange method depends on A and the resolution. 

The error in F x is due to the fact that the derivatives of the velocity field are not considered 
in the boundary conditions. This can be understood by analyzing Eq. (16a). At the steady 
state, and with the approximation that 


fa » + /i” = 



(e„ • V)(e„ • ti) , 


( 21 ) 


Equation (16a) at the top wall becomes 

F x « 2 w 2 p-ze 2 ■ (u bf + u f - 2u w ) , (22) 

c z 

where the substitution of Eq. (10) for / 6 and / 8 has been made. The only term in the 
above equation which has A dependence is u b f. When 0 < A 1/2, F x is independent of 
A, and when 1/2 < A < 1, F x weakly depends on A because u w = 0 in this case [see 
Eqs. (9)]. In the case where F x is obtained by summing over a set of symmetric lattice 
points, cancellations in the summation may further weaken the dependence of F x on A. 
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Table I also shows that for the shear stress based on taking the derivative of the velocity, 
the loss of accuracy is quite significant for small values of A (< 0.05) for r = 0.6. For other 
values of A (> 0.3), the accuracy is comparable with that of F x . However, as shown in 
Fig. 5(a), the accuracy of pv^\ y=H based on the near-wall velocity derivative deteriorates 
as the relaxation time r increases (from 0.51 to 1.6). To see the cause of the increasing error 
in pv^-\ y =H, Fig. 5(b) shows dimensionless wall velocity, u w /u c , obtained by a three-point 
second-order Lagrangian extrapolation of the near wall velocity u x (y ) as a function of r. 
The increasing slip velocity u w on the wall with the increasing relaxation time r was also 
observed in Ref. [16]. It is the result of increasing particle mean free path that causes the 
deviation of the kinetic solution from the hydrodynamic solution. It is clear that the poor 
performance of pv^-\ y=H is associated with the increasing error in the near wall velocity 
profile as r increases. Since the stress tensor Ty can be calculated directly from / Q [see 
Eq. (13)] without the need for directly computing velocity derivatives, the force evaluation 
method based on the evaluation of the velocity gradient in the form of Eq. (12) is not 
recommended. 

B Steady uniform flow over a column of cylinders 

For a uniform flow over a column of circular cylinders of radius r and center-to-center 
distance H (see the left part of Fig. 9 for illustration), symmetry conditions for / a ’s are 
imposed at y = ±.Hj 2. Most of the details of flow field simulation can be found in Ref. [1]. 
The Reynolds number is defined by the diameter of the cylinder d as Re = Ud/ v, where U 
is the uniform velocity in the inlet. It must be noted that for a consistent determination of 
the force, the upstream boundary must be placed far upstream. A shorter distance between 
the cylinder and the boundary will result in higher drag. In this study, it is placed at about 
20 radii to the left of the center of the cylinder. Reducing the distance between boundary 
and the cylinder to 12.5 radii while keeping the rest of the computational parameters fixed 
would increase the drag coefficient by about 1.8% at Re = 100. The downstream boundary 
is located about 25 - 30 radii behind the cylinder to allow sufficient wake development. The 
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simulation is terminated when the following criterion based on the relative 1 , 2 -norm error in 
the fluid region is satisfied, 

\\u(xi,t+ 1) - u{xi,t ) || 2 

g _ g.eo 

a ||u(a:i,i + l )|| 2 

\ Xi€fl 

In this case, e = 10 -6 was chosen for both Re = 10 and 100. 

Following Fornberg [26], the drag coefficient over a circular cylinder of radius r is defined 
as 



pU 2 r' 


(24) 


Figure 6(a) compares Cp obtained from: momentum exchange method, surface stress in- 
tegration, and finite difference result of Fornberg [26] using a vorticity-stream function for- 
mulation at Re = 100, H/r = 20, and radius r ranging from 2.8 to 13.2. For r > 8, both 
methods of momentum exchange and the stress integration give satisfactory results for Cp 
in comparison with the value of 1.248 given in Ref. [26]. This adds credence to the validity 
of Eq. (15) for evaluating the total force on a body. The values of Cp from the momentum 
exchange method have a little less variation than that from the stress integration. Accepting 
an error of less than 5%, the reliable data for Cp can be obtained, using the momentum 
exchange method, for r > 5. That is, 10 lattices cross the diameter of the cylinder are 
necessary to obtain reliable values of the force. This is consistent with the finding by Ladd 
[15]. In the range of 5 < r < 7, the stress integration method gives more scattered result 
than the method of momentum exchange. For smaller radius, i.e., coarser lattice resolution, 
while both methods give poor results (due to insufficient resolution), the stress integration 
yields much larger errors. 

Figure 6(b) compares Cp obtained from the methods of momentum exchange and the 
stress integration for Re = 10. The momentum exchange method seems to gives a converged 
result at larger r (> 8). Based on the data for r > 8, an average values of Cp « 3.356 is 
obtained. In contrast, the stress integration method has a larger scattering than the large 
r result from the momentum exchange method even for r > 8. Averaging over the results 
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for r > 8, the stress integration gives Cd ~ 3.319. The difference between converged 
results of two methods is about 1%. For r less than or around 5, the scattering in C D 
from the stress integration method is much larger than that in the momentum exchange 
method. The conclusions from the comparisons in Fig. 6 are as follows: (i) both methods 
for force evaluation can give accurate results; (ii) the momentum exchange method gives 
more consistent drag; and (iii) in the range of 10 < Re < 100, a resolution of 10 lattices 
across the diameter of the cylinder are needed in order to obtain consistent and reliable drag 
values. In other words, the lattice (grid) Reynolds number Re* should be less than 10 in the 
calculations. 

In the above results presented in Figs. 6(a) and 6(b), the center of the cylinder is placed 
on a lattice grid thus the computational mesh is symmetric with respect to the geometry 
of the cylinder. To test the effect of the mesh symmetry on the accuracy of the force 
evaluation, the calculation of the flow at Re = 10 is repeated with different values of the 
cylinder center offset A x in x direction, or A y in y direction. The radius of the cylinder is 
deliberately chosen to be only 6.4 lattice grids. In order to preserve the mirror symmetry of 
the flow in y-direction, w r e use different boundary conditions for upper and lower boundaries 
(at y ~ ±H). For A x = 0 while varying A y , we use the periodic boundary conditions at 
y = ±H. For A y = 0 while varying A x , we use the symmetric boundary conditions. The 
results of the drag coefficient Cd are presented in Table II. The variation of Cd due to the 
change of the center of cylinder offset from a grid point is less than 1%. We notice that 
the variation in Cd due to A y is larger than that due to A x . This is precisely because of 
the asymmetry of the mesh due to A y / 0, while the offset in rr-direction maintains the 
flow symmetry in y-direction. This asymmetry due to A,, =4 0 results in the change of the 
lift coefficient from O(10 -12 ) to 0(1O -3 ). It is our observation that the accuracy of the 
force evaluation schemes used here is dictated by that of the boundary conditions at the 
solid walls. The error due to the symmetry of the computational mesh with respect to the 
geometry of an object is well bounded. This is also observed in other independent studies 
[27, 28]. 

It is worth noting that the w r all shear stress in the channel flow obtained by using the 
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method of momentum exchange has a relative error proportional to the resolution across 
the channel. For a resolution of 10 - 20 lattices across the diameter considered here, the 
relative error in the drag appears, however, smaller than in the channel flow case. At 
Re = 100, with r > 10, the average value of the drag obtained by using the method of 
momentum method has a 1.7% relative error comparing with Fornberg’s data [26]. If the 
boundary layer thickness is estimated roughly to be 3 x 2r j \/Re ~ 6, there are only about 
6 lattices across the boundary layer over which the velocity profile changes substantially. 
Based on the insight from the channel flow result, it is possible that the deviatoric shear 
stresses on the surface of the cylinder that are effectively incorporated in the method of 
momentum exchange suffer comparable level of error as in the channel flow. The effective 
error cancellation over the entire surface of the body may have contributed to the good 
convergence behavior in the drag shown in Figs. 6(a) and 6(b). 

C Flow over an asymmetrically placed circular cylinder in channel with vortex shed- 
ding 

Schafer and Turek [29] reported a set of benchmark results for a laminar flow over a 
circular cylinder of radius r that is asymmetrically placed inside a channel. In the present 
study, r = 12.8 is used and the center of the cylinder coincides with a grid point. The 
distance from the center of the cylinder to the upper wall and lower wall is h+ = 4.2 r and 
h- = 4.0 r, respectively. This results in A + = 0.76 for the upper wall and A_ = 0.2 for 
the lower wall, respectively. The channel inlet has a parabolic profile and it is placed at 4 
radii upstream of the cylinder center according to the specification of the benchmark test 
[29]. This results in A = 0.2 for the inlet boundary. A zeroth-order extrapolation for f a is 
used at the exit boundary which is located 40 radii downstream of the cylinder center. Thus 
there are a total of 564 x 105 square lattices in the flow field. For Re = 2rU jv — 100 based 
on the average inlet velocity U, the use of relaxation time r = 0.55 requires U = 0.095. 

At this Reynolds number, the flow becomes unsteady and periodic vortex shedding is 
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observed. Figures 7(a) and 7(b) compare the lift coefficient 



and the drag coefficient C D [see Eq. (24)] with the benchmark results in Ref. [29]. We first 
note that the present numerical value of Strouhal number St = 2 r /UT is 0.300, where T is 
the period of the lift curve. This agrees very well with the range of Cl values (0.2995 - 0.305) 
given in Ref. [29]. We note that the difference in C L (t) between the momentum exchange 
method and the surface stress integration method is indiscernible graphically. For the drag 
coefficient CT>(i), it is interesting to note that although there is about 0.25% difference 
between the results given by momentum exchange method and the surface stress integration 
method, both methods of force evaluation give two peaks in the Cd{ t) curves. Physically, 
these two peaks in Co(t) curve correspond to the existence of a weaker vortex and a stronger 
vortex alternately shed behind the cylinder. The difference in the strength of the vortices 
results from the difference: h + /r — 4.2 and h-/r = 4.0 in the passages between the cylinder 
and the channel walls. There is no report on the occurrence of these two peaks in Ref. [29]. 
Instead, a range of the maximum Co (from 3.22 to 3.24) by different researchers was given. 
The present value of the higher peak is well within the range. A further refined computation 
of the present problem using a multi-block procedure [30] with r = 40 in the fine grid region 
yield nearly the same results for Co(t) and Cx(f). 

D Pressure driven flow in a circular pipe 

The steady state flow field was obtained by using D3Q19 model with r = 0.52 [2], Eq. (15) 
is used to evaluate the force on the boundary points along the circumference of the pipe over 
a distance of one lattice in the axial direction. The resulting axial force, F x , is, equivalently, 
the force given by T w (2nr8 x ) where t w is the wall shear stress and r is the pipe radius. For 
a fully developed flow inside a circular pipe, the exact fluid shear stress at the pipe wall is 
given by 

r® xact (27rr) = 7rr 2 ^ . (25) 
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We examine the normalized axial force, 


Figure 8 shows the normalized coefficient r) over a range of r: 3.5 - 23.5. Except for r < 5, 
r] is rather close to 1. It was noticed in Ref. [2] that the accuracy of LBE solution for the 
pipe flow is not as good as that for the two-dimensional channel flow due to the distribution 
of values of A around the pipe. The accuracy of the drag is dictated by the accuracy of the 
flow field if the force evaluation method is exact. For the pipe flow, the error in F x results 
from the inaccuracy in the flow field and the errors in the force evaluation scheme based on 
momentum exchange (as seen in the previous section for the two-dimensional channel flow 
case). For r > 5, the largest error in F x is about 3.5% and it occurs at r = 15.5. Again, there 
is no systematic error in F x . Given the complexity of the boundary in this three-dimensional 
flow, the results shown in Fig. 8 are satisfactory in the sense that it adds further credence 
to the momentum exchange method for force evaluation. 


E Steady uniform flow over a sphere 

To limit the computational effort, a finite domain of —H/2 < y < Hj 2 and —Hf 2 < 
z < H/2, with H/r = 10 is used to compute the flow past a sphere of radius r (see Fig. 9). 
Two cases are considered: (a) the flow past a single sphere, and (b) the flow over an two- 
dimensional array of spheres (all located at x — 0) with the center of the spheres forming 
square lattices. In the former case, the boundary conditions at j y = 1 (y = H/2 corresponds 
to j y — 2) for fa s are given by the following linear extrapolation 

fatix, 1 ,jz) = 2 fa(jx, 2, jz) ~ f Q {j x , 3J Z ) . (27) 

The velocity at j y = 2 is set as 

u(j x ,2J z ) = u(j x ,3,j z ) . (28) 

Similar treatment is applied at y = H/2 and z = ±H/2. In the latter case, symmetry 
conditions are posed on / Q ’s at j y = 1 by using the values of f a 's at j y = 3 (see Ref. [1] for 
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the two-dimensional case). At the inlet, a uniform velocity profile is imposed at j x — 1.5 
(halfway between the first and second lattices). The upstream boundary is located at 7.5 
radii to the left of the sphere center in all simulations. 

For flow over a sphere, the drag coefficient is often expressed as 

c ° = -t-£H = < 29 > 

lp[/ 2 7 rr 2 Re 6nrU pv 

where (j) accounts for the non-Stokesian effect of the drag. For two types of the boundary 
conditions at ( y — ±H/2 and z = ±.£7/2), (j) a denotes the non-Stokesian correction for the 
case where the symmetry conditions are imposed at ( y = ±£7/2 and z = ±i£/2) and (j) ^ 
denotes the results for the case where the extrapolation for f Q is used at (y — ±H/2 and 
z = ±H/ 2) in order to simulate the unbounded flow. 

Figure 10(a) shows the non-Stokesian coefficient (j ) <*, for r = 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 5.1, 
5.2, 5.4, 5.6 and 5.8, for H/r = 10 at Re = 10. The relaxation time is r = 0.7. With this 
range of r, the number of the boundary nodes on the surface of the sphere increases roughly 
by a factor of (5.8/3) 2 « 3.74; the actual counts of the boundary nodes x b gives a ratio 
2370/546 = 4.35. The largest difference is 1.9% between r = 3.0 and r - 3.2 which have 
the least resolution in the cases investigated. For a uniform flow over an unbounded sphere, 
an independent computation using finite difference method based on the vorticity-stream 
function formulation with high resolution gives a drag coefficient 4> ~ 1.7986 at Re = 10. 
The largest difference between this result and the LBE results is 1.36% at r — 3.2. If the 
LBE data for the drag is averaged over the range of r, one obtains <j> ~ 1.8086 which differs 
from 1.7986 by 0.54%. Hence, the LBE solutions based on 3.0 < r < 5.8 give quite consistent 
drag force. Figure 10(b) shows the non-Stokesian correction factor <f> s for a uniform flow 
over a planar array of spheres for 3.0 < r < 5.8 and H/r = 10, at Re = 10. It is important 
to note that with the improvement of the surface resolution by a factor of 4.35, there is little 
systematic variation in </> s (r ). The largest deviation from the average value, ~ 1.963, is 
1.1% at r = 5.0. It is clear that the LBE solution gives reliable fluid dynamic force on a 
sphere at r ~ 3.5 for a moderate value of Re. The set of data for <f> s is inherently more 
consistent than that for <j> oo since the symmetry boundary condition can be exactly specified 
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at y — ±Hj 2 and 2 = ±Hf 2, while the extrapolation conditions given by Eqs. (27) and (28) 
do not guarantee the free stream condition at y = ±Hj 2 and z = ±H/2. "Vet, both (j > 00 and 
( fi s exhibit remarkable self-consistency from coarse to not-so-coarse resolutions. 

IV. CONCLUSIONS 

Two methods for evaluating the fluid force in conjunction with the method of lattice 
Boltzmann equation for solving fluid flows involving curved geometry have been examined. 
The momentum exchange method is very simple to implement. It is shown in the channel 
flow simulation that momentum exchange method is not an exact method. The error in the 
wall shear stress is inversely proportional to the resolution. In two- and three-dimensional 
flows over a bluff body, it can give accurate drag value when there are at least 10 lattices 
across the body at Re ~ 100. The method of integrating the stresses on the surface of the 
body gives similar result when there is sufficient resolution but a much larger uncertainty 
exists when the resolution is limited in comparison with the method of momentum exchange. 
In addition, this method requires considerably more efforts in implementing the extrapo- 
lation and integration on the body surface in comparison with the method of momentum 
exchange. The method of momentum exchange is thus recommended for force evaluation 
on curved bodies. 

This work is supported by NASA Langley Research Center. R. Mei also acknowledges 
partial support of the Engineering Research Center (ERC) for Particle Science and Tech- 
nology at the University of Florida, the National Science Foundation (EEC-9402989), and 
Industrial partners of the ERC. 

APPENDIX A: LBE MODELS IN TWO AND THREE DIMENSIONS 

The nine-velocity (or 9-bit) LBE model on a two-dimensional square lattice, denoted as 
the D2Q9 model, has been widely used for simulations of two-dimensional flows. For three- 
dimensional flows, there are several cubic lattice models, such as the fifteen- velocity (D3Q15), 
nineteen- velocity (D3Q19), and twenty-seven-velocity (D3Q27) models, which have been 
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used in the literature [12]. All these models have a rest particle (with zero velocity) in the 
discretized velocity set {e Q |c>! = 0, 1, . . . , , (b — 1)}. For athermal fluids, the equilibrium 
distributions for the D2Q9, D3Q15, D3Q19, and D3Q27 models are all of the following form 


[8,9] 


yK eq ) 


w a p 


3 ( s 9 , v2 ^2 
1 + -(e Q • w) + — (e Q • u) - 


2c 2 


u 


(Al) 


where w a is a weighting factor and e a is a discrete velocity, c = 6 x /6t is the unit speed, and 
5 X and 8 t are the lattice constant and the time step, respectively. The discrete velocities for 
the D2Q9 models are 

"(0,0), c* = 0 

(±1, 0) c, (0, ±l)c, a = 1, 3, 5, 7 (A2) 

(±1, ±l)c, a = 2, 4, 6, 8, 

and the values of the weighting factor w Q are 


I, « = 0 


w * = < §, a = 1, 3, 5, 7 


[ a = 2, 4, 6, 8. 

For the D3Q19 model, the discrete velocities are 

(0,0), a = 0 

(±1, 0, 0) c, (0, ±1, 0) c, (0, 0, ±l)c, a = 1-6 

[ (±1, ±1, 0) c, (0, ±1, ±l)c, (±1, 0, ±l)c, a = 7-18, 
and the weighting factor w a is given by [9] 


w„ = 


3 ’ 


- ) i 


a = 0 


a = 1-6 


(A3) 


(A4) 


18 ’ 

36 ’ a = ’ 

The discrete velocity sets {e Q } for the D2Q9 and D3Q19 models are shown in Fig. 1. 
The density and velocity can be computed from f a : 


(A5) 



(A6a) 

a q 

(A6b) 
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The speed of sound of the above LBE models is 


1 

C s = —^=C 

s/z 

and the equation of state is that of an ideal gas such that 


p = c]p. 


(AT) 


The viscosity of the fluid is 

v = c)\ 

for the discrete velocity model of Eq. (2). It should be noted that the equilibrium dis- 
tribution function fa ^ is in fact a Taylor series expansion of the Maxwellian /O) [8, 9]. 
This approximation of fjf^ in algebraic form by making the LBE method valid only in the 
incompressible flow limit u/c — f 0. 

Equation (2) is often discretized in space x and time t into the lattice Boltzmann equation 

f a (Xi + e a St,t + 5t) - f a (xi,t) = --[/«(**, t) - fa q) (xi,t)} , (A8) 

r 

where r = A /8 t . For this LBGK model [3, 4], the viscosity in the Navier-Stokes equation 
derived from the above lattice Boltzmann equation is 



The —1/2 correction in the above formula for u comes from the second order derivatives of 
f a when f Q {xi + e a S t ,t + S t ) in Eq. (A8) is expanded in a Taylor series in u. This correc- 
tion in v makes the lattice Boltzmann method formally a second order method for solving 
incompressible flows [9]. Obviously, the physical and computational stabilities require that 
t > 1 / 2 . 
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FIG. 1: Discrete velocity set {e Q }. (top) Two-dimensional nine-velocity (D2Q9) model, 
(bottom) Three-dimensional nineteen-velocity (D3Q19) model. 
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FIG. 2: Layout of the regularly spaced lattices and curved wall boundary. The circles (o), 
discs (•), shaded discs (•), and diamonds (O) denote fluid nodes, boundary locations (*„), 
solid nodes which are also boundary nodes (a: b ) inside solid, and solid nodes, respectively. 



FIG. 3: Distribution of the pressure coefficient C P on the surface of a 2D circular cylinder 
of radius r = 6.6, and center-to-center distance Hjr = 10. The stanation point is located at 
6 = 180°. The result is obtained with r = 0.6 and Re = 40. 
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FIG. 4: The channel flow configuration in the LBE simulations with an arbitrary A. 
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FIG. 5: The LBE simulations of the channel flow, with A = 0.2, 1/3, 0.5, and 0.7. The 
pressure drop is d x p = —1.0 x 10 -6 in lattice units, (a) Ratio between the wall force, 
pvdyUxly^H, evaluated by using Eq. (18), and the exact value — = —Hd x pf 2, given by 
Eq. (19) as a function of r. (b) Normalized wall slip velocity u w /u c as a function of r. 







FIG. 6: The drag coefficient for a uniform flow past a column of cylinders over a range of 
radius r. (a) Re = 100. The dashed line indicate the Cd value of Ref. [26] (Cd = 1.242); 
and (b) Re = 10. The dashed lines indicate the values of Cd averaged over 4 largest radii. 
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FIG. 7: The 2D flow past a cylinder asymmetrically placed in a channel. The variations of 
lift and drag coefficients as function of time t (after an initial run time of t 0 ) are compared 
with the benchmark results in Ref. [29]. The dashed horizontal lines indicate the upper and 
lower bounds in Ref. [29]. The solid and dashed curves are the results obtained by using 
momentum exchange and stress integration, respectively, (a) The lift coefficient Cz,(i). Note 
that the results obtained by using the two methods are indistinguishable on the graph, (b) 
The drag coefficient Co(t)- 
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FIG. 8: The ratio rj between the tangential force F x on the pipe and its exact value 
(nr 2 dp/dx) over a range of pipe radius r. 




FIG. 9: Computational domain for the uniform flow past a sphere of radius r. The dashed 
lines indicate boundaries of computational domain, (left) Unbounded domain in xy plane, 
and (right) bounded domain in yz plane. 
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FIG. 10: Flow past sphere. Variation of the non-Stokesian correction factor 0 = 

—F x /6nrUpv as a function of sphere radius r at Re = 10. The dashed lines are values 
of 4>(r) averaged over r. (a) The flow past a single sphere in an unbounded field (H/ r = oo). 
(b) The flow past a planer array of spheres (H/r = 10). 
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A 

x 10 5 ■ 

-F x x 10 5 

~Txy% X 10 5 

-C^\y=H x 1° 5 

-F 

r y 

P 

0.01 

1.601 

1.6333 

1.6010 

3.5294 

0.3333 

0.3333 

0.02 

1.602 

1.6333 

1.6020 

2.5555 

0.3333 

0.3333 

0.03 

1.603 

1.6333 

1.6030 

2.2309 

0.3333 

0.3333 

0.04 

1.604 

1.6333 

1.6040 

2.0685 

0.3333 

0.3333 

0.05 

1.605 

1.6333 

1.6050 

1.9710 

0.3333 

0.3333 

0.1 

1.610 

1.6333 

1.6100 

1.7760 

0.3333 

0.3333 

0.2 

1.620 

1.6333 

1.6200 

1.6781 

0.3333 

0.3333 

0.25 

1.625 

1.6333 

1.6250 

1.6583 

0.3333 

0.3333 

0.3 

1.630 

1.6333 

1.6300 

1.6451 

0.3333 

0.3333 

0.3333 

1.633 

1.6333 

1.6330 

1.6385 

0.3333 

0.3333 

0.35 

1.635 

1.6333 

1.6350 

1.6357 

0.3333 

0.3333 

0.4 

1.640 

1.6333 

1.6400 

1.6285 

0.3333 

0.3333 

0.5 

1.650 

1.6333 

1.6500 

1.6184 

0.3333 

0.3333 

0.6 

1.660 

1.6333 

1.6600 

1.6214 

0.3333 

0.3333 

0.7 

1.670 

1.6333 

1.6700 

1.6244 

0.3333 

0.3333 

0.8 

1.680 

1.6333 

1.6800 

1.6274 

0.3333 

0.3333 

0.9 

1.690 

1.6333 

1.6900 

1.6305 

0.3333 

0.3333 

0.95 

1.695 

1.6333 

1.6950 

1.6321 

0.3333 

0.3333 

0.99 

1.699 

1.6333 

1.6990 

1.6335 

0.3333 

0.3333 


TABLE I: Comparison of fluid stresses at y = H in a two-dimensional pressure driven 
channel flow with dp/dx = —1.0 x 10 -6 in the lattice units, N y = 35 and r = 0.6 as a 
function of A. Column 2, — given by Eq. (19); Column 3, — F x given by Eq. (16a); 
Column 4, — t"®^ given by Eq. (17); Column 5, —pv^-\ y =H Eq- (18); Column 6, — F y given 
by Eq. (16b); Column 7, pressure p obtained in the simulation. 
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a t 

= 0, periodic BC at y - 

= ±H 


Ay 

0 

0.2 

0.4 

0.6 

0.8 

Co 

3.3780 

3.3755 

3.3576 

3.3636 

3.3755 


Ay ‘ 

= 0, symmetric BC in y 

= ±H 


Ax 

0 

0.2 

0.4 

0.6 

0.8 

C D 

3.3745 

3.3844 

3.3847 

3.3838 

3.3860 


TABLE II: The effect of the symmetry of computational mesh on the force evaluation for the 
steady uniform flow over a column of cylinders. The Reynolds number Re = 10 (t = 0.7), 
and the radius of the cylinder r = 6.45 x . The variation of Co due to the change of the center 
of cylinder offset from a grid point is less than 1%. 
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ABSTRACT 


Compared to Navier-Stokes equation-based approach, 
the method of lattice Boltzmann Equation (LBE) offers 
an alternate treatment for fluid dynamics. However the 
LBE method often employs certain regular lattices to 
maintain a compact and efficient computational 
procedure. This limitation makes it difficult to perform 
flow simulations when there is a need for high 
resolution near the body and/or there is a far-field 
boundary. To resolve these difficulties, a multi-block 
method is developed. In this method, the flow field is 
divided into blocks. In each block, the grid is uniform 
and the grid size is chosen according to the desired 
resolution. In this paper, an accurate interface treatment 
between neighboring blocks is derived to ensure the 
continuity of mass, momentum, and stresses across the 
interface. Several test cases are employed to 
demonstrate that the present multi-block method can 
greatly improve the accuracy and computational 
efficiency. 

1. BACKGROUND OF THE LATTICE 
BOLTZMANN METHOD 


Recently, there has been much progress in 
developing and employing the method of the lattice 
Boltzmann equation (LBE) [1-3] as an alternative, non- 
traditional computational technique for solving 

complex fluid dynamic systems [4-5]. In an NS 
equation-based macroscopic method for computational 
fluid dynamics (CFD), the macroscopic variables of 
interest, such as velocity u and pressure p, are obtained 
by solving the Navier-Stokes (NS) equations [6-8]. In 
the LBE approach, one solves the kinetic equation for 
the particle mass distribution function /*,£,?) and the 
macroscopic quantities (such as mass density p and 
momentum density fm) can then be obtained by 
evaluating the hydrodynamic moments of the 

distribution function/. 

A popular kinetic model is the Boltzmann equation 


** Graduate student * Associate professor, f Professor and Chair, 


with the single relaxation time approximation [9]: 

| -+(-Vf =- 4 (/-/ (0> ) (1 > 

dt A 

where E, is the particle velocity, f 0> is the equilibrium 
distribution function (the Maxwell-Boltzmann 
distribution function), and X is the relaxation time. The 
mass density p and momentum density pu are the 
hydrodynamic moments of the distribution function/. 

To solve for / numerically, Eq. (1) is first 
discretized in the velocity space \ using a finite set of 
velocities {£„} without affecting the conservation laws 
[ 5 , 9 - 11 ], 

%■+&•%— (2) 
dr A 

In the above, fjx, t) = fix, £» t) is the distribution 
function associated with the direction \ a and ^ is 


equilibrium distribution function of the ct-th discrete 
velocity. The 9-bit (or 9-velocity) square lattice model, 
which is also called Q9D2 model (Fig. 1) has been 
widely used for simulating 2-D flows. For Q9D2 
model, we use e a io denote the discrete velocity set and 
we have 
£o = 0 , 

(cos((a - 1)# / 4), sin((a - 1 )n f 4)) 


for a=l, 3, 5, 7, 

g or 1 -v/2(cos((a - l)n / 4), sin((ar - 1)7T / 4)) 


for a=2, 4, 6, 8 (3) 

where c = &x/5 t, 5x and St are the lattice constant and 
the time step size, respectively. The equilibrium 
distributions for Q9D2 model (as well as for some of 
the 3-D lattice models) are in the form of 




i+4-( e „-«) 2 + 

c 


(4) 


where is the weighting factor given by 


4/9, 

a= 0 

1/9, 

a = 1,3,5, 7 

1/36, 

or = 2, 4, 6, 8. 


With the discretized velocity space, the density and 
momentum flux can be evaluated as 


P-ifa-ifF' 

k - 0 *=0 
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and 


(6a) 
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pu = £ e a f a = ie a f% q) (6b) 

*=l 

The speed of sound of this model is c s =cl>j3 and the 
equation of state is that of an ideal gas, 

p = pc] ( 7 ) 

Equation (2) can be further discretized in space and 
time. The completely discretized form of Eq. (1), with 
the time step 5 1 and space step eJ5t y is: 
f a (Xi + ejk, t + St) r)= 

( 8 ) 

T 

where T =A/8t, and is a point in the discretized 
physical space. The above equation is the lattice 
Boltzmann equation [1-3] with Bhatnagar-Gross-Krook 
(BGK) approximation [9]. The viscosity in the NS 
equation derived from Eq. (8) is 

v = (2t-\)c s & ( 9 ) 

This choice for the viscosity makes formally the LBGK 
scheme a second order method for solving 

incompressible flows [10,11]. The positivity of the 
viscosity requires that x> V4. Equation (8) can be solved 
as: 

collision step: 

/„(*,. t) = / B (*i.O ■ — [ /« ( Jf / . t) - fa^ ( x i • t)] dOa) 

T 

streaming step: 

fM+e a Stj + St) = Ux^) dOb) 

where ~ denotes the post-collision state of the 
distribution function. It is noted that the collision step 
is completely local and the streaming step takes very 
little computational effort. Equation (10) is explicit, 
easy to implement, and straightforward to parallelize. 

Those inherent advantages of the LBE method 
require the use of regular lattice (such as a square lattice 
or hexagonal lattice) and that the lattice spacing be 
equally distributed. This is in direct contrast to the 
many finite difference/volume/element methods in 
which body fitted coordinates can be used and the grid 
stretching can be easily applied. However, it should be 
noted that there has also been growing interest in the 
macroscopic method to employ the Cartesian grid for 
complex flow problems [12,13]. A challenge of such 
Cartesian grid approaches is to offer high resolution 
near the body and to place the outer boundary far away 
from the body. In order to use the regularly spaced 
lattice while developing the capability to place the outer 
boundary far away, it is desirable to divide the 
computational domain into a number of blocks within 
which a fixed, constant lattice spacing can be used and 
Eq. (10) is implemented in each block as in the standard 
LBE method. Again, such an approach has been 
actively employed in the macroscopic methods with 
both Cartesian and curvilinear coordinates. 


This paper describes a multi-block strategy for the 
LBE method. In each block, constant value of 5x=5y=8t 
is used. The information exchange on the interface 
between the neighboring blocks of different lattice 
spacing 5x for the primary variables fds is implemented 
to ensure the mass conservation and the continuity of 
stresses between blocks. A lid-driven cavity flow is 
computed using a single block with uniform grid and 
the present multi-block method. The results are 
compared with published benchmark results . A channel 
flow with a parabolic velocity profile at the inlet over 
an asymmetrically placed cylinder at Re=100 (based on 
the average incoming velocity) is computed next using 
the multi-block method. Finally, flow over NACA0012 
airfoil at Re=500-5000 is computed. The present study 
shows that the multi-block strategy can greatly improve 
the computational efficiency of the LBE method. 


2. BASICS OF THE MULTI-BLOCK 
STRATEGY IN THE LBE METHOD 


To illustrate the basic idea, a two-block system (a 
coarse and a fine, as shown in Fig. 2) is considered in 
the derivation for the interfacial information exchange. 
The ratio of the lattice space between the two-grid 
system is 

m - 5x c / Sxf (11) 

For a given lattice size 5x, the viscosity of the fluid is 

v = (2r-l)8 x c/ 6 (12) 

In order to keep viscosity v, and thus Re , the same in 
the entire flow field involing different lattice sizes, the 
relation of relaxation times, ry on the fine grid and x c on 
the coarse gird, must obey the following relation: 

Tf = — + m(r c “) (13) 

J 2 2 

for c=l. To keep the variables and their derivatives 
continuous on interface between two systems of 
different grids, consistent, accurate relationship 
between the two grid systems must be developed. 

The Chapman-Enskog expansion gives, 

/ a (r,t) =f ( a q) (x,t)+ r) + ... (14) 




nf(*9) Df (eq) 

= -A £lSL- s - 

Dt Dt 

It is noted that 

/-(!)_ f f(«9) _ Anon-eq) 

Ja ~ Ja ~Ja ~ Ja 

is the non-equilibrium part of the distribution function 
based on which the deviatoric stresses are evaluated. 

The collision step (Eq. (10a)) gives 


(15) 

(16) 
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/«(*,-. 0= a--)fa(x i ,t) + -f^ eq \x i ,t) (17) 

T T 

Sustituting (9) into (17) leads to 

T 

+ti' ) (X i , t) + ...}+- f^\x„ t) 

T 

= f^ eq) (x i ,t)+— /®(*,,t) + ... (18) 

T 

Denoting the coarse-grid results with superscript c and 
fine-grid result with superscript/, the post-collision step 
gives 


?c _ t(eq,c) ^ /r * + ... 


_ f\eqs) 

J a ~ J a T 


(19) 


Similarly, 

7/ = /i^ /) +^-l/i , ’ /) +- (2°) 

r / 

Since the velocity and density must be continuous on 
the interface between the two grids and fa^ from 
equation (4), it is seen that 

fa qX) = fa qJ) (2D 

To maintain the continuity in the deviatoric stresses, in 
the 2-D case, 

r ij = ( 1 ~h')i\.fa-fa“ ,) 1(«. 

a - 1 


‘■ca’-aj ^ ^ a® \j ) 


2r 

= u~) I -\e a e a Sy) (22) 

2r 2 

it is obvious that one requires 

(1 _JL) f{non-eq,c) _ \ ) j>(non-eqJ) (23) 


2r t 


It i 


or 




(!/') - 


( 1 - 


2r 


-> fS J) 


or 


(If) 


T f 

Substituting (24) into (19) one obtains 

fa = fa qX)+mlL ^fa’ f) +■■■ 

T f 

Using Eqs. (20, 21), the above becomes 


f c = f(eqj) 
J a Ja 


+ m 


r c -\ r 


f 


= f Kt 

J a 


( eqj ) 


+ m- 


r /~ 

r / _1 


[// 

• Iff -f ( a qJ) ) 


(24) 


(25) 


(26) 


In transferring the data from the coarse grids to the fine 
grids, one similarly obtains 


J f - 

Ja 


fS 


{eq,c) 


jLlzL 

m(T c - 1) 


[fa 


■/r- c) ] 


(27) 


On the interface between two blocks, there are m 
values of f& needed for each fa q ' c ^ and /<£ . Thus, 
spatial and temporal interpolation procedures for the 
values of fi eq,c) and fa on the fme -g rid lattice is 


used to complete the evaluation of . There are two 
kinds of interpolation are used along the interface: 
interpolation in space and interpolation in time. To 
eliminate the possibility of asymmetry caused by 
interpolations, a symmetric, 4-point spline fitting is 
used for spatial interpolation. We found that it is very 
important to maintain the symmetry in the interpolation 
along the interface. For example if a 3-point Lagrangian 
interpolation is used for spatial interpolation, the 
asymmetry caused by the interpolation can generate a 
lift coefficient of 0.01 for steady flow over cylinder at 
Re=20. And if flow is inherently unsteady, this 
asymmetry could exaggerate the extent of the 
unsteadiness over long time. For temporal interpolation, 
a 3-point Lagrangian interpolation is used. 


3. RESULTS AND DISCUSSIONS 

In all three cases considered here, the boundary 
condition for f a in the solid region near a wall is 
obtained using the formulations given in Ref. [14] for a 
curved geometry. 

3.1 Lid-driven cavity flow 


The lid-drive cavity flow has been extensively used 
as a benchmark solution to test the accuracy of a 
numerical method in spite of the fact that two singular 
points at the upper corner require high resolution to 
obtain smooth stresses near the corner points. To assess 
the LBE results, the bench mark solution of Ref. [15] 
are used for comparison. 

The computations are carried out using a single- 
block with uniform lattice (129x129) with the walls 
placed halfway between lattices and a multi-block 
whose layout is shown in Fig. 3. Near the two upper 
corner points, the grid resolution is increased by factor 
4. The relaxation time is t c = 0.56 for coarse-grid block 
and T f =0.74 for the fine-grid block. The upper wall 
velocity is U=0.0155945. The initial condition for the 
density is unity and that for velocity is zero. The 
streamlines shown in Fig. 4 are obtained from the single 
block solution and the pattern is not discernable from 
those of the multi-block solution. The positions of the 
centers of the primary vortices are (0.6154, 0.7391) and 
(0.6172, 0.7390) for uniform grid and multi-block 
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respectively, comparing with the value (0.6172, 0.7344) 
from reference [15]. The u- and v-components of the 
velocity along the vertical line and horizontal line 
through geometry center are shown in Fig 5a and 5b, 
respectively, along with the bench mark result. It is 
seen that while the single block method with 129x129 
lattices can capture most of physical variables 
accurately, the multi-block method gives more accurate 
velocity profiles. Fig. 5 shows the pressure contour 
from the single-block computation. Because of the 
singularity at the upper corners, the density contour has 
a very large region of spatial oscillations due to the 
insufficient resolution near the singularity. Fig. 7 shows 
the pressure contours obtained from the multi-block 
solution. Significant improvement in the smoothness of 
the solution for the pressure field over that of the single 
block solution is observed. If the finest grid size is used 
in the entire region, it will require a significantly longer 
CPU time and much larger computer memory. 

In an NS solver for incompressible flows, because 
the decoupling of thermodynamic pressure and velocity 
field, it is crucial to maintain the mass conservation of 
the entire flow domain. This issue becomes more 
critical when the multi-block method is used [8, 16]. 
Also for incompressible flows, the pressure is arbitrary 
up to a constant. Hence coupling the pressure term 
while maintaining the mass flux conservation is very 
important. Generally speaking, it is difficult to maintain 
simultaneously the continuity of mass, momentum, and 
stresses across the interface between neighboring 
blocks because interpolations are applied separately 
along the interface to evaluate mass flux, momentum 
flux, and pressure. Hence it is often difficult to use 
large grid size ratio (m) between two neighboring 
blocks. In multi-block LBE method, the continuities of 
mass and stresses are ensured through the use of 
Equations (26-27). The most important point is that 
interpolations are only applied to /' s along the interface 
and this automatically ensures the consistency in the 
transfer of various flux terms across the interface. 

To validate the above arguments, pressure, shear 
stress, mass flux and momentum flux near the block 
interfaces are examined next. Fig. 8 shows the local, 
enlarged view of the pressure contour around an 
interface corner point indicated by the circle in Fig. 7. 
Clearly, the pressure is rather smooth across the 
interface with the coarse-to-fine grid size ratio of m- 4. 
Figs. 9-11 show the contours of shear stress, mass flux, 
and momentum flux pu x 2 . It is seen that these physical 
quantities are all smooth across the interface. 

To demonstrate this issue more clearly, 
macroscopic physical quantities on one part of the 
interface (i.e. line A-B in Fig. 3) are plotted in Figs. 12- 
17. After streaming step there is no physical value on 
the interface for the fine grid. Here we use second order 
extrapolation to get fine grid value on interface. Figs. 


12-15 show that mass and momentum flux match very 
well between fine- and coarse-grid. Fig. 16 shows the 
shear stress. In most part of the interface two sets of 
values agree very well with each other. The discrepancy 
appears near the upper wall. It is noted that for in the 
fine-grid blocks, the top moving wall is located half- 
way between two horizontal, fine-grid lattices with a 
distance of AjSxf =0.5 Sx f . In the coarse-grid block, the 
distance between the wall to the nearest lattice in the 
fluid region is A c &c c =0.5 Sx f =0.5 Sx c /4=0.125 Sx c for 
m= 4. This mismatch (A f *A C ) will result in different 
errors in the boundary condition for ff s. This 
subsequently affects the accuracy of the shear stress 
near the corner of block and wall. The same problem 
also appears in Fig. 17 for pressure. 

3.2 Channel flow over an asymmetrical placed 
cylinder at Re=100 

Schafer & Turek [17] reported some benchmark 
results for a laminar flow over a circular cylinder 
placed asymmetrically inside a channel. The cylinder 
has a radius of 0.1m and is asymmetrically placed in the 
channel. The center to the upper wall distance is 
hjr= 0.21m and the center to the lower wall distance is 
hJn=> 0.20m. In the LBE computation, r=5 lattice is used 
in the coarse grid system and the coarse-to-fine lattice 
spacing ratio is m= 4. The coarse-grid block has a total 
of 220x42 lattices. The relaxation times are t c = 0.52 
for and Tf =0.58. The channel inlet has a parabolic 
velocity and is located about 4 radii upstream of the 
cylinder center. A zeroth-order extrapolation for/ a is 
used. The Reynolds number based on the average inlet 
velocity and the diameter of the cylinder is Re^ 100. 

At this Reynolds number, the flow becomes 
unsteady and periodic vortex shedding is observed. The 
numerical value of Strouhal number is 0.300 and it 
agrees very well with the value (0.2995-0.305) in 
reference [17]. An instantaneous streamline plot is 
shown in Fig. 18 after the dynamically periodic solution 
is established. The drag and lift are shown in Fig. 19. 
The unsteady characteristics of the flow agree well with 
the reported results in [17] and comparable with a 
similar computation using lattice Boltzmann method 
[18]. 

3.3 Steady flow over NACA0012 airfoil 

The NACA 0012 airfoil (Fig. 20) is a popular wing 
model, which has been used extensively. Flow fields at 
Re=500, 1000, 2000, and 5000 are presently computed 
with the multi-block LBE scheme. Fig. 21 shows the 
entire computational domain and the schematic of the 
multi-block arrangement. There are 150 lattices (grids) 
along the chord in the finest block. At the inlet, upper, 
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and lower boundaries, the equilibrium boundary 
condition is used. At the downstream boundary a zeroth 
order extrapolation for/’s is used. 

Fig. 22 shows the density contour, streamlines and 
velocity vector of the converged solution at Re=2000 
and zero angle of attack. To investigate the effect of 
grid resolution, two sets of grid systems are used for the 
flow field at Re=500: a fine grid system and a coarse 
grid system (with resolution reduced by a factor of 2 in 
every block from the fine grid system). Fig. 23 shows 
the velocity profiles at (x-Xle)/L=0.06 where L is the 
chord length and x LE is the location of the leading edge. 
The two sets of velocity profiles agree well with each 
other, although the fine grid solution appears to have 
smoother u-component velocity profile, as expected. 

Fig. 24 compares the drag coefficient C d between 
the present LBE simulation and those calculated from 
Xfoil [19]. It can be seen that two sets of results agree 
with each other very well for the entire range of 
Reynolds numbers investigated in this study. 

It is also noted that at /te=500, the present value of 
Q=0.1761 compare very well with the results reported 
in Ref. [20]: Q=0.1762 obtained using a Navier-Stokes 
equation-based finite difference method and Q=0. 1717 
using Powerflow code developed by EXA Corporation, 
which is based on the lattice Boltzmann equation 
method. In addition, the present simulation for the 
symmetrical flow at Re=500 gives a lift coefficient of 
I C L \ <6 X 10' 14 . Ref. [20] reported C L =1.1 5x1 O' 7 using 
an NS equation-solver and C L =2.27xl0‘ 4 using EXA’s 
Powerflow code. This suggests that the present multi- 
block code preserve the symmetry very well. 

Finally, it is worth pointing out that there is a 
tremendous saving in the computational cost using the 
multi-block method in LBE simulations. There are three 
different sizes of grids used for the NACA0012 airfoil 
simulation. There are 513x65=33345 fine grids, 23300 
intermediate grids with m= 4, and 37407 coarse grids 
(with m=8 in reference to the finest grids). This gives a 
total of about 9.4x1 0 4 grids in the entire domain. If the 
fine grid system is used in the entire domain, the 
number of the grids would be N x xN y ^= 2849x577- 
1.64X10 5 which is 18 times more than in the multi- 
block case. This represents a saving of 18 times in the 
memory. Furthermore, since 8f=5jc=8y in the LBE 
simulation, one time step marching in the coarsest grid 
system (m=8) requires 2 sweeps in the intermediate grid 
blocks and 8 sweeps in the finest grid blocks. The ratio 
of the computational efforts required to carry out a 
single-block simulations to that for a multi-block 
simulation for a given period of physical time would be 
1.64xl0 5 x8/(33345x8+23300x2+37407) - 38. Clearly, 
more saving can be achieved if more blocks of different 
sizes are used. 


4. CONCLUDING REMARKS 

A multi-block strategy is developed for the lattice 
Boltzmann method. The interface condition is derived 
to ensure the mass conservation and stress continuity 
between neighboring blocks. Favorable computational 
results are obtained in three test cases. There is a 
significant potential for the multi-block strategy in the 
LBE method in the aerodynamics application since the 
boundary conditions at infinity and on the wall can both 
be reconciled. 
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Fig. 2 Interface characteristics between two 
blocks of different lattice spacing. 



Fig 4 Streamlines in the cavity flow at 
Re=100. 
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Fig 5a Comparison of u- velocity along the 
vertical line through geometric center. 



Fig 6 Pressure contours in the cavity flow 
from the single-block LBE simulation. 




360 380 400 


X 

Fig. 8 Enlarged view of pressure contour in 
the circled region in Fig. 7 near the 
intersection of three blocks. The figure 
demonstrates that the block interface and 
corner are well handled. 



Fig 5b Comparison of v-velocity along the 
horizontal line through geometric center. 



Fig.7 Pressure contours in the cavity from 
multi-block LBE solution. (For the circled 
region, see Fig. 8) 



Fig. 9 Shear stress contour. Solid and dash 
lines represent positive and negtive values, 
respectrively. 
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Fig. 10 Contour of x-component mass flux pu x . 
Solid and dash lines represent positive and negtive 
values, respectively.. 


Fig. 1 1 Contour of momentum flux in the x 
2 

direction p-u . 


firm grid 

■o— - coarse gird 


— • fine grid 

— c - - coarse gird 


Fig. 12 The x-component of the mass flux pu x I{fbU) Fig. 13 The x-component of the mass flux 

on the interface AB. In Figs 12-15, fo - 1 and pu i(poU) on the interface AB. 

U=0.0155945. 


coarse grid 
fine grid 



Fig. 14 The x-component of the momentum flux , 
2 2 

pu x l p 0 U , on the interface AB. 


Fig. 15 The y-component of the momentum flux, 
2 

pUxUylPoU , on the interface AB. 
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Fig. 16 Shear stress f{fxU ! H ) on the interface Fig. 17 Pressure on the interface AB. 
AB. 



Fig. 18 Instantaneous streamlines for channel flow over an asymetrically placed cylinder at 
Re=100. 




Fig. 19 Unsteady drag and lift coefficients on the cylinder 
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Fig. 20 NACA 0012 airfoil. 
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Fig. 21 Block and lattice layout for flow over NACA 0012. The lattice 
spacing is reduced by a factor 16 for graphical clarity. 



550 600 650 700 750 800 850 900 050 1000 

X 

Fig. 22 Streamlines, pressure contour, velocity vector for flows over 


NACA 0012 airfoil at Re=2000. 





AIAA 2000-2614 



Fig. 23 Grid-independence test of the velocity 
profiles near the leading edge at (x-x le )/L=0.06 for 
flow over NACA0012 airfoil at Re=500. 



Fig. 24 Comparison of C d between the present 
simulation and Xfoil calculation as a function of Re 
for flow over NACA0012 airfoil. The straight line is 
the slope according to the laminar boundary layer 
theory. 
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The generalized hydrodynamics (the wave vector dependence of the transport coefficients) of a generalized 
lattice Boltzmann equation (LBE) is studied in detail. The generalized lattice Boltzmann equation is con- 
structed in moment space rather than in discrete velocity space. The generalized hydrodynamics of the model 
is obtained by solving the dispersion equation of the linearized LBE either analytically by using perturbation 
technique or numerically. The proposed LBE model has a maximum number of adjustable parameters for the 
given set of discrete velocities. Generalized hydrodynamics characterizes dispersion, dissipation (hyperviscosi- 
ties), anisotropy, and lack of Galilean invariance of the model, and can be applied to select the values of the 
adjustable parameters that optimize the properties of the model. The proposed generalized hydrodynamic 
analysis also provides some insights into stability and proper initial conditions for LBE simulations. The 
stability properties of some two-dimensional LBE models are analyzed and compared with each other in the 
parameter space of the mean streaming velocity and the viscous relaxation time. The procedure described in 
this work can be applied to analyze other LBE models. As examples, LBE models with various interpolation 
schemes are analyzed. Numerical results on shear flow with an initially discontinuous velocity profile (shock) 
with or without a constant streaming velocity are shown to demonstrate the dispersion effects in the LBE 
model; the results compare favorably with our theoretical analysis. We also show that whereas linear analysis 
of the LBE evolution operator is equivalent to Chapman -Enskog analysis in the long-wavelength limit (wave 
vector k=0), it can also provide results for large values of k. Such results are important for the stability and 
other hydrodynamic properties of the LBE method and cannot be obtained through Chapman-Enskog analysis. 


PACS number(s): 47.10.Tg, 47.1 l.+j, 05.20.Dd 

I. INTRODUCTION 

The method of lattice Boltzmann equation (LBE) is an 
innovative numerical method based on kinetic theory to 
simulate various hydrodynamic systems [1-3]. Although the 
LBE method was developed only a decade ago, it has at- 
tracted significant attention recently [4,5], especially in the 
area of complex fluids including multiphase fluids [6-11], 
suspensions in fluid [12], and viscoelastic fluids [13,14], The 
lattice Boltzmann equation was introduced to overcome 
some serious deficiencies of its historic predecessor: the lat- 
tice gas automata (LGA) [15-17]. The lattice Boltzmann 
equation circumvents two major shortcomings of the lattice 
gas automata: intrinsic noise and limited values of transport 
coefficients, both due to the Boolean nature of the LGA 
method. However, despite the notable success of the LBE 
method in simulating laminar [18-21] and turbulent [22] 
flows, understanding of some important theoretical aspects 
of the LBE method, such as the stability of the LBE method, 
is still lacking. It was only very recently that the formal 
connections between the lattice Boltzmann equation and the 
continuous Boltzmann equation [23-25] and other kinetic 
schemes [26] were established. 
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In this work we intend to study two important aspects of 
the LBE method which have not been systematically studied 
yet: (a) the dispersion effects due to the presence of a lattice 
space and (b) conditions for stability. We first construct a 
LBE model in moment space based upon the generalized 
lattice Boltzmann equation due to d’Humieres [27]. The pro- 
posed model has a maximum number of adjustable param- 
eters allowed by the freedom provided by a given discrete 
velocity set. These adjustable parameters are used to opti- 
mize the properties of the model through a systematic analy- 
sis of the generalized hydrodynamics of the model. General- 
ized hydrodynamics characterizes dispersion, dissipation 
(hyperviscosities), anisotropy, lack of Galilean invariance, 
and instability of the LBE models in general. The proposed 
generalized hydrodynamic analysis enables us to improve the 
properties of the models in general. The analysis also pro- 
vides us better insights into the conditions under which the 
LBE method is applicable and comparable to conventional 
computational fluid dynamics techniques. 

Furthermore, from a theoretical perspective, we would 
like to argue that our approach can circumvent the Chapman- 
Enskog analysis to obtain the macroscopic equations from 
the LBE models [27,13,14]. The essence of our argument is 
that the validity of the Chapman-Enskog analysis is entirely 
based upon the fact that there are two disparate spatial scales 
in real fluids: the kinetic (mean-free-path) and the hydrody- 
namic scales the ratio of which is the Knudsen number. 
When the LBE method is used to simulate hydrodynamic 
motion over a few lattice spacings, there is no such separa- 
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tion of the two scales. Therefore, the applicability of 
Chapman-Enskog analysis to the LBE models might become 
dubious. Under the circumstances, analyzing the generalized 
hydrodynamics of the model becomes not only appropriate 
but also necessary. 

It should also be pointed out that there exists previous 
work on the generalized hydrodynamics of the LGA models 
[28-32] and the LBE models [33]. However, the previous 
work only provides analysis on nonhydrodynamic behavior 
of the models at finite wavelength, without addressing im- 
portant issues such as the instability of the LBE method or 
providing insights into constructing better models. In the 
present work, by using a model with as many adjustable 
parameters as possible, we analyze the generalized hydrody- 
namics of the model so that we can identify the causes of 
certain nonhydrodynamic behavior, such as anisotropy, and 
lack of Galilean invariance, and instability. Therefore, the 
analysis shows how to improve the model in a systematic 
and coherent fashion. 

This paper is organized as follows. Section II gives a brief 
introduction to the two-dimensional (20) nine-velocity LBE 
model in discrete velocity space. Section III discusses the 
generalized LBE model in moment space. Section IV derives 
the linearized lattice Boltzmann equation from the general- 
ized LBE model. Section V analyzes the hydrodynamic 
modes of the linearized evolution operator of the generalized 
LBE model, and the generalized hydrodynamics of the 
model. The dispersion, dissipation, isotropy, and Galilean 
invariance of the model are discussed. The eigenvalue prob- 
lem of the linearized evolution operator is solved analytically 
and numerically. Section VI analyzes the stability of the 
LBE model with Bhatnagar-Gross-Krook (BGK) approxima- 
tion, and compares with the stability of the LBE model pre- 


sented in this paper. Section VII discusses the correct initial 
conditions in the LBE simulations, and presents numerical 
tests of shear flows with discontinuities in the initial velocity 
profile. Section VIII provides a summary and concludes the 
paper. Two appendices provide additional analysis for varia- 
tions of the LBE models. Appendix A analyzes a model with 
coupling between density p and velocity u, and Appendix B 
analyzes the LBE models with various interpolation 
schemes. 

II. 21) NINE-VELOCITY LBE MODEL 

The guiding principle of the LBE models is to construct a 
dynamical system on a simple lattice of high symmetry 
(mostly square in 2D and cubic in 3D) involving a number of 
quantities that can be interpreted as the single-particle distri- 
bution functions of fictitious particles on the links of the 
lattice. These quantities then evolve in a discrete time ac- 
cording to certain rules that are chosen to attain some desir- 
able macroscopic behavior that emerges at scales large rela- 
tive to the lattice spacing. One possible “desirable 
behavior” is that of a compressible thermal or athermal vis- 
cous fluid. (To simplify the analysis, in this work we shall 
restrict our analysis to the athermal case.) We shall demon- 
strate that the LBE models can satisfactorily mimic the fluid 
behavior to an extent that the models are indeed useful in 
simulating flows according to the similarity principle of fluid 
mechanics. For the sake of simplicity, we confine our discus- 
sions here to two-dimensional space. The extension to three- 
dimensional space is straightforward, albeit tedious. 

A particular two-dimensional LBE model considered in 
this work is the nine-velocity model. In this model, space is 
discretized into a square lattice, and there are nine discrete 
velocities given by 


e 


a 


(0,0), tt = 0, 

(cos[(a— 1 )7r/2],sin[(a- l)7r/2])c, a= 1-4, 

(cos[(2a — 9) 7r/4],sin[(2a — 9)7t/4]) V^c, ot= 5—8, 


(0 


where c = S x IS t is the unit of velocity and S x and S t are the 
lattice constant of the lattice space and the unit of time (time 
step), respectively. From here on we shall use the units of 
S x = 1 and S t - 1 such that all the relevant quantities are di- 
mensionless. The above discrete velocities correspond to the 
particle motion from a lattice node r ; to either itself, one of 
the four nearest neighbors (a=l-4), or one of the four 
next-nearest neighbors (a=5-8). This model can easily be 
extended to include more discrete velocities and in space of 
higher dimensions, thereby including further distant neigh- 
bors to which the particles move in one time step. Neverthe- 
less, “hopping” to a neighbor on the lattice induces inherent 
limitations in the discretization of velocity space. 

For the particular model discussed here, nine real num- 
bers describe the medium at each node r } of a square lattice: 


The number f a can be considered as the distribution function 
of velocity e„ at location r y (and at a particular time r). The 
set {/„} can be represented by a vector in R 9 that defines the 
state of the medium at each lattice node: 

l/(r 7 )) = (/o./i,...,/g) T . (2) 

Once the vector |/(r ; )) is given at a point r y in space, the 
state of the medium at this point is fully specified. 

The evolution of the medium occurs at discrete times t 
= nS t (with <5,*= 1). The evolution consists of two steps: (1) 
motion to the relevant neighbors (modeling of advection) (2) 
redistribution of the {f a } at each node (modeling of colli- 
sions). These steps are described by the equation 


{/ a (r,)|a=0,l,...,8}. 


/a( r j + e a' /+ 1 ) = fa( r j ’ 0 + ^a(/)- 


(3) 
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The above equation is the so-called lattice Boltzmann equation (LBE). The lattice Boltzmann equation can be rewritten 
concise vector form: 

|/(r;+e,,f+l)H/(rj')> + |A/>, 


in a 
(4) 


where the following notations are adopted: 

|/(r j +e a ,r+l)}-[/o(r ; + e 0 ,r + l),/ ) (r ; +e 1 ,r+l), •• ,/ 8 (r ; +e 8 ,t+l)] T , 

|A/>=[fto(/),n,(/) fW)] T . 


(5a) 

(5b) 


so that |/(r J + e 0 ,/+ 1)) is the vector of a state after advec- 
tion, and |A f) is the vector of the changes in |/) due to 
collision fl. 

The advection is straightforward in the LBE models. The 
collisions represented by the operator fl may be rather com- 
plicated. However, fl must satisfy conservation laws and be 
compatible with the symmetry of the model (the underlying 
lattice space). This might simplify fit considerably. One 
simple collision model is the BGK model [34,2,3]: 

fta— < 6) 

where t is the relaxation time in units of time step 8 t (which 
is set to be 1 here), and / a eq) is the equilibrium distribution 
function that satisfies the following conservation conditions 
for an athermal medium: 


p=SA eq) =2/ a . 

a a 

(7a) 

i=2 ejr = 2 ej a . 

(7b) 


a a 


where p and u are the (mass) density and the velocity of the 
medium at each lattice node, respectively. For the so-called 
nine-velocity BGK model, the equilibrium is usually taken as 


, aP 


l+3(e„-u)+-(e a -u) 2 -^u 2 


( 8 ) 


where w 0 = 4/9, h'i i2 . 3 , 4 = 1/9, and w 5 , 6 , 7 , 8 = 1/36. 

Some shortcomings of the BGK model are apparent. For 
instance, because the model relies on a single relaxation pa- 
rameter r, the Prandtl number must be unity when the model 
is applied to thermal fluids, among other things. One way to 
overcome these shortcomings of the BGK LBE model [2,3] 
is to use a generalized LBE model which nevertheless retains 
the simplicity and computational efficiency of the BGK LBE 
model. 


III. MOMENT REPRESENTATION 
AND GENERALIZED 2D LBE 

Given a set of b discrete velocities, {eja = 0,l, - . . >(b 
-1)} with corresponding distribution functions, {f a \a 
= 0,1, . . . ,(fc — 1)}, one can construct a fo-dimensional vec- 
tor space R 6 based upon the discrete velocity set, and this is 
usually the space mostly used in the preceding discussion of 
the LBE models. One can also construct a space based upon 
the (velocity) moments of {/„}. Obviously, there are b inde- 
pendent moments for the discrete velocity set. The reason in 


favor of using the moment representation is somewhat obvi- 
ous. It is well understood in the context of kinetic theory that 
various physical processes in fluids, such as viscous trans- 
port, can be approximately described by coupling or interac- 
tion among “modes” (of the collision operator), and these 
modes are directly related to the moments (e.g., the hydro- 
dynamic modes are linear combinations of mass, and mo- 
menta moments). Thus the moment representation provides a 
convenient and effective means by which to incorporate the 
physics into the LBE models. Because the physical signifi- 
cance of the moments is obvious (hydrodynamic quantities 
and their fluxes, etc.), the relaxation parameters of the mo- 
ments are directly related to the various transport coeffi- 
cients. This mechanism allows us to control each mode in- 
dependently. This also overcomes some obvious deficiencies 
of the usual BGK LBE model, such as a fixed Prandtl num- 
ber, which is due to a single relaxation parameter of the 
model. 

For the nine-velocity LBE model, we choose the follow- 
ing moments to represent the model: 

|p) = (l,l,U,U,l,U) T . (9a) 

|e) = ( — 4,— 1- L- 1 - 1,2,2,2,2) t , (9b) 

|e) = (4,2,2, 2,2,1, 1,1,1 ) T , (9c) 

1;,) = (0,1,0,- 1,0,1,— 1,- 1,1 ) T , (9d) 

| </x ) = (0,- 2,0,2, 0,1,- 1,- 1,1 ) T , (9e) 

|y v ) = (0,0, 1,0,— 1,1,1, 1, 1) T » (90 

|, 7 ,) = (0,0,- 2,0,2, 1,1- 1,-1 ) T , (9g) 

|p„) = (0,l,- LI - 1,0,0,0,0) T , (9h) 

\ Pxy ) = (0,0,0, 0,0,1,- 1,1,- 1) T . (90 

The above vectors are represented in the space V=R 9 
spanned by the discrete velocities {e a }, and they are mutu- 
ally orthogonal to each other. These vectors are not normal- 
ized; this makes the algebraic expressions involving these 
vectors which follow simpler. Note that the above vectors 
have an explicit physical significance related to the moments 
of {/„} in discrete velocity space: |p) is the density mode; 
|e) is the energy mode; |e) is related to energy square; \j x ) 
and | j y ) correspond to the x and y components of momentum 
(mass flux); | q x ) and Ip^) correspond to the x and y compo- 
nents of energy flux; and | p xx ) and \p xy ) correspond to the 
diagonal and off-diagonal component of the stress tensor. 
The components of these vectors in discrete velocity space 
V=R 9 are constructed as follows: 

|p)a=|e.l°=l. ( 10a ) 
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\e) a = -4|eJ°+3(eL +e oo>)’ 

(10b) 

= 4|e a |°— y {e 2 a _ x + e\ y )+ (e 2 a<x + e 2 a J 2 , 

(10c) 

l/jc) « e a,x ’ 

(10d) 

\q x ) a = [-5\e a \°+3(ei x + el y )]e aa . 

(10e) 

\jy)a~ e a,y > 

(100 

\q y ) a = [-5\e a \° + l(e 2 a . x + el.y)]e a .y’ 

(10g) 

|Pxx)a — e a,x e a,y * 

(10h) 

\Pxy} a~~ e a, x e a, y • 

(10i) 

p=(pI/>=(/1p>* 

(11a) 

e = (e\f) = (f\e). 

(lib) 

e = (e |/> = (/|s). 

(11c) 

j X =(j X \f) = (f\jx)> 

(lid) 


<7,= <<?,!/>= (/IgA ( lle > 

jy = (jy\f) = (f\jy)’ < llf > 

g,=(^l/)=(/kv)- < n 8) 

pxx=(pxxl/>=(/lpxx>, ( llh ) 

P X y=(P*y\f) = (f\P*y)- < 1U > 

Similar to {/„}, the above set of moments can also be con- 
cisely represented by a vector: 

| Q)-={p,e,e,j x ,q x ,j y ,q y >p xx >Px>) T - ( 12 ) 

There obviously exists a transformation matrix M between 
|(?) and |/) such that 

le>-M|/>, (13a) 

|/) = M _1 |(?). (13b) 


In other words, the matrix M transforms a vector in the vec- 
tor space V spanned by the discrete velocities into a vector in 
the vector space M—R^ spanned by the moments of {/< r } ■ 
The transformation matrix M is explicitly given by 


!<Pl \ 


1 1 

1 

1 

1 

1 

1 

1 

1 

l 1 

<4 1 


1 -4 

-1 

- 1 

-1 

-1 

2 

2 

2 

2 

( e l 


4 

-2 

-2 

— 2 

-2 

1 

I 

1 

1 

Oxl 


0 

1 

0 


0 

1 

-1 

- I 

1 

<9x1 

= 

0 

-2 

0 

2 

0 

1 

-1 

-1 

1 

0/1 


0 

0 

1 

0 

-1 

1 

1 

-1 

-1 

(<7yl 


0 

i 

0 

-2 

0 

2 

1 

1 

- 1 

-1 

i (Pxxl 1 


0 

1 

-I 

1 

-1 

0 

0 

0 

0 

\ (Pxyl 1 

I 

\ 0 

0 

0 

0 

0 

1 

-1 

1 

-1 


= (|p)’l e )>l e )’l./x)'kx)'l./.y)’l‘7y)>|Pxx)’IPx>’) ' 


(14) 


The rows of the transformation matrix M are organized in the 
order of the corresponding tensor, rather than in the order of 
the corresponding moment. The first three rows of M corre- 
spond to p, e, and e, which are scalars or zeroth-order ten- 
sors, and they are zeroth-order, second-order, and fourth- 
order moments of {f a }, respectively. The next four rows 
correspond to j x , q x , j y , and q y , which are vectors or first- 
order tensors, and j x and j y are the first-order moments, 
whereas q x and q y are the third-order ones. The last two 
rows represent the stress tensor, which are second-order mo- 
ments and second-order tensors. Again, this can easily be 
generalized to models using a larger discrete velocity set, 
and thus higher-order moments, and in three-dimensional 
space. The main difficulty when using the LBE method to 
simulate a real isotropic fluid is how to systematically elimi- 
nate as much as possible the effects due to the symmetry of 
the underlying lattice. We shall proceed to analyze some 
simple (but nontrivial) hydrodynamic situations, and to make 
the flows as independent of the lattice symmetry as possible. 


I 

Because the medium simulated by the model is athermal, 
the only conserved quantities in the system are density p and 
linear momentum j~ (j x J y )- Collisions do not change the 
conserved quantities. Therefore, in the moment space M, col- 
lisions have no effect on these three quantities. We should 
stress that the conservation of energy is not considered here 
because the model is constructed to simulate an athermal 
medium. Moreover, we find that the nine-velocity model is 
inadequate to simulate a thermal medium because it cannot 
have an isotropic Fourier law for the diffusion of heat. Al- 
though the conserved moments are not affected by collisions, 
the nonconserved moments are affected by collisions, which 
in turn cause changes in the gradients or fluxes of the con- 
served moments, which are higher-order moments. In what 
follows, the modeling of the changes of the nonconserved 
moments is described. 

Inspired by the kinetic theory for Maxwell molecules 
[35], we assume that the nonconserved moments relax lin- 
early towards their equilibrium values that are functions of 
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the conserved quantities. The relaxation equations for the 
nonconserved moments are prescribed as follows: 


e* = e — s 2 [e~ e (eq) ]. 

(15a) 

i — i 
O' 

'“co 

1 

CO 

1 

CO 

II 

* 

CO 

(15b) 

<7? = < 

(15c) 

q*=q y -s ?[<?). -<7> eq) ]’ 

(15d) 

Pxx=Pxx- s i\.Pxx~P { xx\ 

(15e) 

Pxy=Pxy -S&Pxy ~ > 

(15f) 

where the quantities with and without superscript * are post- 
collision and pre-collision values, respectively. The equilib- 
rium values of the nonconserved moments in the above equa- 
tions can be chosen at will provided that the symmetry of the 
problem is respected. We choose 

e(eq) = ^ e ^[a 2 (p\p)p+y 2 ((j x \j x )j x + (j > \j y )j y )] 

1 1 , , 

= J«2P+g y2(j x + Jy)’ 

(16a) 

B (eq) = ^ E ^[ay(p\p)p+ y i ((j x \j x )j x + (j y \jy)jy)] 

1 1 , 

= 4«3 P + ^yiU'x+jy)’ 

(16b) 

(eq) _ (Jx\Jx) 1 

q * {q x \qx) Ux 2 lJ *' 

(16c) 

(ea) (Jy\jy) . 1 

(M Clh 2 Clh ' 

(16d) 


p [? )= y 1 (p J^ (Ox\jx)j 2 x-(jy\jy)jl)= 2 

(16e) 


to solve the mathematically difficult problem to create an 
interparticle collision mechanism for the fictitious particles 
in the LBE models that would give the same eigenmodes of 
the collision operator in the continuous Boltzmann equation. 
However, what can be accomplished is that by carefully 
crafting a simple model with certain degrees of freedom, we 
can optimize large-scale properties of the model in the sense 
that generalized hydrodynamic effects (deviations from hy- 
drodynamics) are minimized. 

The values of the unknown parameters, Cj, a 2t 3 , and 
7 i, 2 , 3 , 4 > shall be determined by a study of the modes of the 
linearized collision operator with a periodic lattice of size 
N x XN y . 

It should be noted that in Eq. (16) the density p does not 
appear in the terms quadratic in j. This implies that the den- 
sity fluctuation is decoupled from the momentum equation, 
similar to an incompressible LBE model with a modified 
equilibrium distribution function [36]: 



P + Po 


3(e„-u)+-(e Q u) 2 --u 2 


(17) 


where the mean density p 0 I s usually set to be 1. The model 
corresponding to the equilibrium distribution function of Eq. 
( 8 ) shall be analyzed in Appendix A. 


IV. LINEARIZED LBE 

We consider the particular situation where the state of the 
medium is a flow specified by uniform and steady density p 
(usually p=l, so the uniform density may not appear in 
subsequent expressions) and velocity in Cartesian coordi- 
nates \=(V x ,V y )y with a small fluctuation superimposed: 

|/) = |/ (0, > + |tf>, 08 ) 

where |/* 0) ) represents the uniform equilibrium state speci- 
fied by uniform and steady density p and velocity V 
= (V X >V y ), and |<5/) is the fluctuation. The linearized Bolt- 
zmann equation is 

| Sf( r, + e, ,t- + 1 )) = | Sf( tj ,t ) ) + n<°> | Sf( r, ,t )) (19) 


tea) y l(jx\Jx)(jy\Jy) , . . N \ 

The values of the coefficients in the above equilibria (c ^ , 
<* 2 , 3 * and 71 , 23 , 4 ) will be determined in the next section and 
summarized in Sec. V E. The choices of the above equilibria 
are made based upon inspection of the corresponding mo- 
ments given by Eqs. (10), or the physical significance of 
these moments. Note that in principle q x and q y can include 
terms involving third-order terms in terms of moment, such 
as j\ and j x p xx [14], and e can include fourth-order terms. 
Nevertheless, for the nine-velocity model, these terms of 
higher order are not considered because either they do not 
affect the hydrodynamics of the model significantly, or they 
lead to some highly anisotropic behavior which is undesir- 
able in the LBE modeling of hydrodynamics. 

Clearly, LBE modeling of fluids is rather different from 
real molecular dynamics. Therefore, it is not necessary to try 


where H (0) is the linearized collision operator: 




l/>= L/< 0> > 




( 20 ) 


In the moment space M, the linearized collision operator can 
be easily obtained by using Eqs. (15) and (16): 


&,3 a 

where Q a and j Q a ) y <* = 0,1, . . . ,(fe- 1) are the moments 
defined by Eqs. (11) and the corresponding vectors in V 
-R 9 defined by Eqs. (9); is the change of the moment 
due to collision given by Eqs. (15); |f?) = |£? (0> ) * s the vector 
of all moments at the uniform equilibrium state [see Eq. (12) 
for the definition of | £)]. Obviously the linearized collision 




(Qa\e*) JQfi 


(21) 


le)=le (0) > 
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operator C depends on the uniform state specified by density p and velocity V-(V x ,V y ), upon which the perturbation \Sf) is 
superimposed. Specifically, for the nine-velocity model. 


C= 


I 0 

s 2 a 2 /4 

530-3/4 

0 

0 

0 

0 

0 

0 


0 

“*2 

0 

0 

0 

0 

0 

0 

0 


0 

0 

“*3 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 \ 



*272^/3 

0 

5 272 V3 

0 

0 

0 1 



•5374^/3 

0 

5374V/3 

0 

0 

0 



0 

0 

0 

0 

0 

0 



s$c | / 2 

^5 

0 

0 

0 

0 


(22) 

0 

0 

0 

0 

0 

0 



0 

0 

5 7 Cj / 2 


0 

0 



358 71V, 

0 

-35 8 y,y y 

0 

-5g 

0 

j 



35973 V y /2 

0 

35 9 73V,/2 

0 

0 

-S 9 I 



to I Sf) is 

! 

|5/(k,f+I))» 

L|5/(k.r)>. 

(28) 


| Sq), and | Sq) = M| Sf ) . The change of the perturbation due 
to collisions is linearly approximated by |Ae) = C|<5f?) in 
the moment space M spanned by {|g a )| a = 0,1, . . . ,(£> 
- 1)}. This change of state in discrete velocity space V is 
\Af) = M~ ] C\Sg). Therefore, Eq. (19) becomes 


I Sf(Tj+e a ,/+ 1)) = | Sf(Tj ,t )) + M -1 CM| Sf(rj ,t)). 


(23) 


where 

L=A _l [l + M -1 CM], 
is the linearized evolution operator. 

V. MODES OF LINEARIZED LBE 


(29) 


In Fourier space, the above equation becomes 

A|<?/(k,f+ 1 )) = [! + M _1 CM]|<5/(k,r))» (24) 

where A is the advection operator represented by the follow- 
ing diagonal matrix in discrete velocity space V=R 9 : 

A afi - exp(ie„-k)<$ a/3 , (25) 

where 8 a/3 is the Kronecker delta. It should be noted that for 
a mode of wave number k=(k x ,k y ) in Cartesian coordi- 
nates, the advection operator A in the above equation can be 
written as follows: 

A=diag(l,p,q, \lp,\.lq,pq,qlpAlpq,plq)> (26) 


A. Hydrodynamic modes and transport coefficients 

The evolution equation (23) is a difference equation that 
has a general solution: 

\G(rj ,t=l))= K™K"z l \G 0 ), (30) 

where m and n are indices for space (r,- = mx+ny), and x 
and y are units vectors along the x axis and y axis, respec- 
tively; |G 0 ) is the initial state. We can consider the particular 
case of a periodic system such that the spatial dependence of 
the above general solution can be chosen as 

|<5/)= exp(-/k-r ; + zr)|G(r ; ,r)). (31) 


where 

p = e lk q = e lk y. (27) 

The advection can be decomposed into two parts, along two 
orthogonal directions, such as x axis and y axis in Cartesian 
coordinates: 

A(k x ) = A(k x ,k y =0) = diag{l,p,lMp,l,p,lfpA/p,p), 

A(k y ) = A(k x =0,k y ) = d\ag{ \,\,q,\,\lq,q,qAlqAlq), 

and AC/tjt) and A (k y ) commute with each other: 

A= A{k x )A(k y ) = A{k y )A(k x ), 

i.e., the advection operation can be applied along the x direc- 
tion first, and then along the y direction, or vice versa. The 
linearized evolution equation (24) can be further written in a 
concise form: 


By substituting Eqs. (30) and (31) into the linearized LBE 
(28), we obtain the following equation: 

z|G 0 ) = L|Go), (32) 

The above equation leads to the dispersion relation between 
z and k: 

det[L-zl] = 0, (33) 

which determines the transport behaviors of various modes 
depending on the wave vector k. The solution of the above 
eigenvalue problem of the linearized evolution operator L 
provides not only the dispersion relation, but also the solu- 
tion of the initial value problem of Eq. (28): 

b 

|<S/(k,r)) = L'|<5/(k,0))= 2 z‘ a \z a )(z a \8f (k,0)), 

a= 1 
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where \z a ) is the eigenvector of L with eigenvalues z a in 
discrete velocity space V. 

The eigenvalue problem of Eq. (33) cannot be solved ana- 
lytically in general, except for some very special cases. Nev- 
ertheless, it can be easily solved numerically using various 
packages for linear algebra, such as lapack. For small k, it 
can be solved by a series expansions in k. The only part of L 
that has k dependence is the advection operator A. Therefore, 
we can expand A" 1 in L: 

K= A" 1 = K c0) + K (1) (k) + K (2) (k 2 ) + ■ ■ • + K {n \k n ) + 

(34) 

where K (,l) depends on k n : 

K (-»'*■ *■>"*«/»• (35) 

When k=0, the eigenvalue problem of the (bXb) matrix 
L (0) = (1 + M -1 CM) can be solved analytically. There exists 
an eigenvalue of 1 with threefold degeneracy, which corre- 
sponds to three hydrodynamic (conserved) modes in the sys- 
tem: one transverse (shear) and two longitudinal (sound) 
modes. It is interesting to note that when k=(7r,0) or k 
= (0,77-), L also has an eigenvalue of - 1 , which corresponds 
to the checkerboard mode, i.e., it is a conserved mode of L 2 . 
Being a neutral mode as far as stability is concerned, it will 
be necessary to study how it is affected by a mean velocity 
V. Thus we shall have to analyze the model for k ranging 
from 0 to 77, which the standard Chapman-Enskog analysis 
cannot do. 

The hydrodynamic modes at k=0 are 

\q t )= COS 8\j x ) — sin 6\j y )~ \j T ) , (36a) 

le±> = |p>-(cos 6\ j x )+ sin 0\j y ))=\p)±\j L )' (36b) 

where 6 is the polar angle of wave vector k. For finite k, the 
behavior of these hydrodynamic modes depends upon k. In 
two-dimensional space, these linearized hydrodynamic 
modes behave as follows [37]: 

l£?r(0) = Zt\Qt(®)) 

= exp[ — ik(g V cos <p)t ] 

Xexp( — vk 2 t)\Q T (0 )), (37a) 

le±(0>=*±le±(o» 

= exp [±ik(c s ±gV cos cf>)t] 

X exp[ - ( v!2 + £)k 2 t] \ Q ± (0)), (37b) 

where v and £ are the shear and bulk viscosity, respectively; 
the coefficient g indicates whether the system is Galilean 
invariant (that g=l implies Galilean invariance); c s is the 
sound speed; V is the magnitude of the uniform streaming 
velocity of the system V=(V X ,U > ); and (p is angle between 
the streaming velocity V and the wave vector k. The 
Galilean-coefficient g( k) is similar to the g factor in the 
Frish-Hasslacher-Pomeau lattice gas automata [15-17], 
which also determines the Galilean invariance of the system. 


The transport coefficients and the Galilean-coefficient are re- 
lated to the eigenvalues of L as the following: 

v(k)= -Re[lnz r (k)], (38a) 

k~ 

g(k)V cos <£= — -Im[lnz r (k)], (38b) 

^i/(k) + £(k)=- — Re[lnz ± (k)], (38c) 

^ k L 

c s (k)±g(k)V cos <p~ +-Im[lnz ± (k)], (38d) 

where z r (k) and z ± (k) are the eigenvalues corresponding to 
the hydrodynamic modes of the linearized evolution operator 
L. Since the transport coefficients can be obtained through a 
perturbation analysis, we shall use the following series ex- 
pansion in k: 

v(k)= j/ 0 - V]fc 2 + ... +(-l ) n v n k 2n +---, (39a) 

f(k) = £o-fi k 2 + l) n ( n k 2n +---, (39b) 

c s (k) = C 0 -C l k 2 + . .. +(-l) n C,,* 2 "+---, (39c) 

£(k)=,?o-gi* 2 + ■••+(- D n g n k 2n +---. (39d) 

It should be noted that, in the usual Chapman-Enskog analy- 
sis of LBE models, one only obtains the values of the trans- 
port coefficients at k=0. As we shall demonstrate later, 
higher-order corrections to the transport coefficients (i.e., hy- 
perviscosities) are important to the LBE hydrodynamics, es- 
pecially for spatial scales of a few lattice spacings. 

One possible method by which to solve the dispersion 
relation det[L-zl] = 0 is to apply the Gaussian elimination 
technique using l/s a as small parameters for the noncon- 
served modes (the kinetic modes). Starting from a 9X9 ( b 
X b in general) determinant, we obtain a 3 X 3 determinant 
for the three conserved modes. The elements of this new 
determinant are computed as a series of l/j a and k with the 
necessary numbers of terms to achieve a given accuracy 
when computing the roots of the dispersion equation. 

It should be mentioned that the value of the present tech- 
nique is that it provides a very simple means by which to 
analyze models with various streaming and collision rules 
with as many adjustable parameters as possible to be deter- 
mined later when trying to satisfy either the stability criteria 
or physical requirements to model various hydrodynamic 
systems. Free parameters are the equilibrium coefficients in 
Eqs. (16): c { , a t , and y, ; and relaxation rates . 

B. Case with no streaming velocity ( V— 0) 

We first consider the case in which the streaming velocity 
V=0. To the first order in k, we obtain two solutions of 
Im(ln z ± )= +ikc s with 
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These are the sound modes supported by the medium. At the 
next order, we obtain modes with Re(ln z T )= ~ . To en- 

force isotropy we need to have 


1 1/1 l\(c,+4) 

^ 2 _2 U 2/ (2-c,) ’ 

such that the 0 dependence in v 0 vanishes, 

(2-c,) / 1 1\ 


(41) 


(42) 


which can be interpreted as the shear viscosity of the me- 
dium in the limit k= 0 (measured in basic units of space and 
time). For the sound modes, we also find an attenuation rate 

Re(lnz ± )=-(V2+£ 0 )* 2 where (V 2 + fo) is the longitu- 
dinal kinematic viscosity in a two-dimensional system. The 
bulk viscosity of the model at long-wave-length limit k=0 is 


£o= 


(c, + 10— 12c 2 ) j 1 1 

24 U _ 2 


(43) 


The positivity of the transport coefficients leads to the 
bounds on the adjustable parameters: 


C s = V cos 4>± \Jc] + V 2 cos 2 <£, (48) 

where Vcos (j)=\-k, and k is the unit vector parallel to k. 
This clearly shows that the system obeys Galilean invariance 
only up to first order in V. One way to correct this defect is 
to allow for compressibility effects in the equilibrium prop- 
erties, as shown in Appendix A. The dispersion of sound can 
be computed either analytically, by carrying out the pertur- 
bation expansion in k, or numerically, by solving the eigen- 
value problem for any value of k. The dispersion of sound is 
important when studying the nonlinear acoustic properties of 
the medium. 

Second, the attenuation of the transverse wave depends 
not only on V but also on the direction of the wave vector k. 
In order to eliminate the anisotropy in the V dependence of 
the shear wave attenuation, we must choose 

c, = — 2. (49) 

With the above choice of c, t the shear viscosity in the limit 
of k=0 is given by 

«'oK*2(2-*8)[c j 2 + (1-3 c, 2 )V 2 cos 2 <£] + 3[2(i 8 -5 2 ) 


— 16< a 2 , 


-4<c,<2, 


(44a) + s%{s 2 ~- 2 ) cos2 <£]^ 4 cos2 <f>V[ 6 s 2 s% 


(44b) 


X(V 2 cos 2 <£ + c^)]. 


(50) 


and the bounds on the following relaxation parameters: 

0<5 2 <2, (45a) 


Similarly, from the attenuation of acoustic waves, one ob- 
tains the bulk viscosity (in the limit of k=0) that has a 
complicated dependence on the streaming velocity V: 


0<s g <2. (45b) 

The bounds for a 2 and c { will be further narrowed in the 
following analysis. Based upon the above results of v 0 , £ 0 , 
and c s , it is clear that the model is isotropic at rest (i.e., the 
streaming velocity V=0) and in the limit of k=0. The 
Galilean-coefficient g cannot be determined when the 
streaming velocity V=0. Therefore, the case of a finite 
streaming velocity V is considered next. 

C. Case with a constant streaming velocity V 

As indicated by Eqs. (38), to the first order in k, the three 
hydrodynamic roots of the dispersion equation (z r and z±) 
give the phase g Vcos <£ and the sound speed c s . In order to 
cause the root of the transverse mode (z r ) to have a correct 
phase corresponding to the streaming velocity V, as expected 
for a model satisfying Galilean invariance, i.e., we 

must set 

2 

ri = r3=j- ( 46 ) 

If we further set 

y 2 = 18, (47) 

then we obtain the roots of the sound modes (z T _) which lead 
to the sound speed 


C 0 = (Vcos <K'V 2 cos 2 ^ + c J i {12V 2 [(5 2 -i 8 ) 

+ i2 (s 8 -2)cos 2 0] + (2s 2 — 3 s 2 s 8 + 4j 8 )( 1 -3c 2 )} 

+ 3V 4 cos 2 <£[cos 2 <£(2j 8 -F 8s 2 )"f 6(* y 2"~ 5 8)] 

+ 2 V 2 cos 2 <fr[ 6 (s 2 Ss-S 2 - s%)c 2 s +s s (2- s 2 )] 

+ c 2 s [ 6 V 2 {s 2 -s s ) + s s (2-s 2 ) 
X(2~3c 2 )])/{12j 2 j b (V 2 cos 2 <j> +c 2 )}. (51) 

It is obvious that the streaming velocity V has a second order 
effect on v 0 , and a first-order effect on £ 0 . A careful inspec- 
tion of the above result of £ 0 indicates that the first-order 
effect of V on f 0 can be eliminated by setting c 2 s - 1/3 (or, 
equivalently, a 2 = —8). Furthermore, the second-order effect 
of V on the sound speed and the longitudinal attenuation can 
also be eliminated by using a slightly more complicated 
model with thirteen velocities, as noted by a previous work 
[38]. 

In summary, although all the transport coefficients are 
isotropic in the limit k=0, some undesirable features of the 
LBE models can be clearly observed at the second order in k 
when the streaming velocity V has a finite magnitude. First, 
the acoustic wave propagation is not Galilean invariant. Sec- 
ond, both the shear and the bulk viscosities depend on V. 
Nevertheless, these effects are of second order in V, and can 
be improved to higher order in both k and V by incorporat- 
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ing compressibility into the equilibrium properties of the mo- 
ments (see Appendix A) or using models with a larger ve- 
locity set. 


D. Third-order result 


The analysis in the preceding subsections shows that isot- 
ropy for the hydrodynamic modes of the dispersion equation 
can be attained to the first and second orders in k by care- 
fully adjusting the parameters in the model. In the situation 
with a uniform streaming velocity V parallel to k, we find 
that the third-order term in k for the shear mode is aniso- 
tropic, i.e., 


8 1 = 


2 

34 



2 

— + 
*8 


4i 


V 2 cos 2 4> 


1 

1 

l 

2 \ 

/ , , l\ 



+ — 

1 

cos 4 0— cos~ 6+ ^ 

3 

*8 

s 5\ s 8 1 

11 3/ 


(52) 


The anisotropic term in g j (depending on cos ff) can be 
eliminated if we choose 


5 5 = 3 


(2 s g ) 


(53) 


relaxation parameters s ( in the model as opposed to one in 
the LBE BGK model. Two of them, 5 2 and s % , determine the 
bulk and the shear viscosities, respectively. Also, because 
C[ = — 2, 5 9 = 5 8 [see Eq. (41)]. The remaining three relax- 
ation parameters, 5 3 , s 5 , and j 7 , can be adjusted without 
having any effect on the transport coefficients in the order of 
k 2 . However, they do have effects in higher-order terms. 
Therefore, one can keep values of these three relaxation pa- 
rameters only slightly larger than 1 (no severe over- 
relaxation effects are produced by these modes) such that the 
corresponding kinetic modes are well separated from those 
modes more directly affecting hydrodynamic transport. 

It is interesting to note that the present model degenerates 
to the BGK LBE model [2,3] if we use a single relaxation 
parameter for all the modes, i.e., 5 *= 1 /r, and choose 

a 3 = 4 , (54a) 

y 4 = -18. ( 54b ) 

Therefore, in the BGK LBE model, all the modes relax with 
exactly the same relaxation parameter so there is no separa- 
tion in time scales among the kinetic modes. This may se- 
verely affect the dynamics and the stability of the system, 
due to the coupling among these modes. 


As indicated by Eq. (42), parameter 5 8 is usually chosen 
close to 2 from below in order to obtain a small shear vis- 
cosity (and, consequently, a large Reynolds number). There- 
fore, the preceding expression yields a small value for 5 5 . 
This would lead to an undesirable consequence: Mode | q x ) 
relaxed with the relaxation parameter 5 5 would become a 
quasiconserved mode leading to some sort of viscoelastic 
effect [14]. Therefore, we usually choose to have large 5 5 
such that the advection coefficient of transverse waves has an 
angular dependence for nonzero k in third order in k. That is, 
the physical conservation laws are preserved at the expense 
of the isotropy of the dispersion in third order (and all higher 
orders) in k. 

It should be noted that the value of g has effects on the 
Reynolds number because the time t needs to be rescaled as 

<?'• 

E. Optimization of the model and connection 
to the BGK LBE model 

Among seven adjustable parameters (c l5 or, , and y,) in 
the equilibrium values of the moments in the model [see Eqs. 
(16)], so far only five of these parameters have been fixed by 
enforcing the model to satisfy certain basic physics as shown 
in the preceding analysis: Cj = — 2, a 2 =- 8 , y 1 = y 3 =2/3, 
and 72 = 18. These parameter values are the optimal choice 
in the sense that they yield the desirable properties (isotropy, 
Galilean invariance, etc.) to the highest order possible in 
wave vector k. It should be stressed that the constraints im- 
posed by isotropy and Galilean invariance are beyond the 
conservation constraints — models with only conservation 
constraints would not necessarily be isotropic and Galilean 
invariant in general, as observed in some newly proposed 
LBE models for nonideal gases [39,40,9]. Two other param- 
eters, ar 3 and y 4 , remain adjustable. In addition, there are six 


VI. LOCAL STABILITY ANALYSIS 

The stability of the LBE method has not been well under- 
stood, although there exists some preliminary work [41,42]. 
However, previous work does not provide much theoretical 
insight into either the causes or the remedies for the instabil- 
ity of the LBE method. In the following analysis, a system- 
atic procedure that identifies some causes of instability is 
discussed and illustrated by some examples. 

Our stability analysis relies on the eigenvalue problem for 
the linearized evolution operator L, the dispersion equation. 
For large values of k, one could in principle analyze the 
dispersion equation to higher order by perturbation expan- 
sion. In practice, it is more efficient to compute the roots of 
the dispersion equation numerically. We shall try to identify 
the conditions under which one of the modes becomes un- 
stable: instability occurs when Re(lnz a )<0. 

We have noticed some interesting qualitative properties of 
the dispersion for the nine-velocity model when wave vector 
k is parallel to certain special directions with respect to the 
lattice line. These properties are listed in Table I. These 
qualitative behaviors of the dispersion equation already dem- 
onstrate the strong anisotropy of the dispersion relations dic- 
tated by the lattice symmetry. 

To exhibit the complex behavior of the dispersion equa- 
tion, we compute the roots of the dispersion equation with a 
given set of parameters. Figures 1(a) and 1(b) show the real 
and imaginary parts of the logarithm of the eigenvalues as 
functions of k, respectively. Figure 1 clearly exhibits the 
coalescence and branching of the roots. This suggests a com- 
plicated interplay between the modes of collision operator 
affecting the stability of the model. The asymmetric feature 
of these curves is due to the presence of a constant stream- 
ing. 

The growth rate of a mode \z a ), Re(ln z a ), depends on all 
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k 

Dispersion equation 

Conditions 

( 0 , 0 ) 

[z-l ] 3 = 0 
[z-( 1 -j 2 )] = 0 
[z-(l-r 3 )] = 0 
[z-(l-*s )] 2 = 0 
[z-O -* 8 )] 2 =0 

57 = 55 

(± 1,0)77 
or 

(0,± 1)77 

[z+l] = 0 

[z + (l -s 5 )] = 0 or [z+(l-i 7 )] = 0 
[z + (l-i 8 )] = 0 or [z+(l-i 9 )] = 0 
[z 2 -55 5 z + s 5 - 1] = 0 

+ ?{S2(- s 8 -4 ' s 3)~ 6i 3 i 8 + 9(' s 2 +s 3 + i 8 _2 )} 22 

+ j(r 8 — l)(r 2 (r 3 - 2 ) + s 3 )z 
+ (1— J 2 )(l — i 3 )(l — ^8)] = 0 


(±l,±l)rr 

[z-(1-* 8 )] 2 =0 

[z 2 -5r 5 z + r 5 _1 ] 2 = 0 
[z 3 + 9 (1 li 2 - 3s 3 -9)z 2 
+ £{3(4s 3 — 3)-s 2 (s 3 + 2)}z 
+ (l-s 2 )(l-r 3 )] = 0 



the adjustable parameters: the relaxation parameters, the 
streaming velocity V, and the wave vector k. To illustrate 
this dependence, we consider the BGK LBE model with 
l/r= 1.99. Figure 2 shows the growth rate for the most un- 
stable mode as a function of streaming velocity V and wave 
vector k. For each V, we let k be parallel to V. with a polar 
angle 6 with respect to the x axis. Then we search for the 
most unstable mode in the interval O^k^ir. For the nine- 
velocity BGK LBE model, the unstable mode starts to appear 
above V«=0.07. Figure 2 shows the strong anisotropy of the 
unstable mode: the growth rate significantly depends on the 
direction of k, and the critical value of k at which the un- 
stable mode starts to appear is also strongly anisotropic. We 
also compute the growth rate for the most unstable mode 
with V perpendicular to k, and find that the stability of the 
model is generally qualitatively the same as when V is par- 
allel to k, but is slightly more stable. Generally, we find that 
the transverse mode is more stable than longitudinal modes. 
In many instances we have observed that sound waves 
propagating in the direction of the mean flow velocity V can 
be quite unstable. This instability may be reduced by making 
the first-order V-dependent term in the attenuation of the 
sound waves [£ 0 in Eq. (51)] equal to 0 by choosing c; 
= 1/3, as indicated in the preceding section. It should be 
noted that when the growth rate is infinitesimal, it takes an 
extremely long time for the instability to develop in simula- 
tions. Because the unstable modes we have observed have a 
large wave vector k (small spatial scale), as a practical means 
of reducing the effect of instabilities in LBE simulations, 
some kind of spatial or temporal filtering technique may be 
used in the LBE schemes to reduce small-scale fluctuations 
and thus to limit the development of instabilities. 

It should be pointed out that we do not discuss here the 
influence of boundary conditions that may completely 
change the stability behavior of the model through either 



0 rr/2 TT 

k 



FIG. 1. Logarithmic eigenvalues of the nine- velocity model. 
The values of the parameters are a 2 — — 8, 0:3 = 4, C] — -2, y 1 
= 73 = 2/3, y 2 = 18, and y 4 =-18. The relaxation parameters are 
3 2 =1.64, s 3 = 1.54, j 5 = j 7 =1.9, and s 8 = s 9 =1.99. The streaming 
velocity V is parallel to k with V = 0.2, and k is along the x axis, (a) 
Re(lnz„) and (b) Im(ln z a ). 
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FIG. 2. Growth rate of the most unstable mode for the BGK 
LBE model -In z a vs the streaming velocity magnitude V. The 
relaxation parameter j 8 = l/r= 1.99. The wave vector k is set par- 
allel to the streaming velocity V. For each value of V with a polar 
angle 6 with respect to the x axis, the growth rate is computed in the 
interval 0<k^ ir in k space. Each curve corresponds to the growth 
rate of the most unstable mode with a given V, and k parallel to V 
with the polar angle 6 with respect to the x axis. 

large-scale genuine hydrodynamic behavior or local excita- 
tion of Knudsen modes. 

As previously indicated, the adjustable parameters in our 
model can be used to alter the properties of the model. The 
stability of the BGK LBE model and our model is compared 
in Fig. 3. In this case we choose the adjustable parameters in 
our model to be the same as the BGK LBE model, but main- 
tain the freedom of different modes to relax with different 
relaxation parameters s a . Figure 2 shows that for each given 
value of V, there exists a maximum value of — 1 /r (which 
determines the shear viscosity) below which there is no un- 
stable mode. The values of other relaxation parameters used 
in our model are 5 2 = L63, 5 3 = 1.14, 55 = 57 = 1.92, and 5 9 
= 5 8 = 1/r. Figure 3 clearly shows that our model is more 
stable than the BGK LBE model in the interval 1 . 9 ^ 5 8 



FIG. 3. Stability of the generalized LBE model vs the BGK 
LBE model in the parameter space of V and 5 8 - 1/r. The lines with 
symbols □ and X are results for the BGK LBE model and the 
model proposed in this work, respectively. The region under each 
curve is the stable region in the parameter space of V and 5 8 = 1/r. 
Note that the stability of the BGK LBE model starts to deteriorate 
after 5 8 ^ 1 . 92 , whereas the stability of the proposed generalized 
LBE model remains virtually intact. 


= l/r^l.99. Therefore, we can conclude that by carefully 
separating the kinetic modes with different relaxation rates, 
we can indeed improve the stability of the LBE model sig- 
nificantly. 

VII. NUMERICAL SIMULATIONS 
OF SHEAR FLOW DECAY 

To illustrate the dispersion effects on the shear viscosity 
in hydrodynamic simulations using the LBE method, we 
conduct a series of numerical simulations of the shear flow 
decay with different initial velocity profiles. The numerical 
implementation of the model is discussed next. 

A. Numerical implementation and initial conditions 

The evolution of the model still consists of two steps, 
advection and collision. The advection is executed in discrete 
velocity space, namely, to |/(x, 0 )» but not to the moments 
|p(x,f)). However, the collision is executed in moment 
space. Therefore, the evolution involves transformation be- 
tween discrete velocity space V and moment space M, simi- 
lar to Fourier transform in the spectral or Galerkin methods. 
The evolution equation of the model is 

|/(x+ e„ <5, + <5,)) = |/(x,0) + M' 1 S[ I e (x,0) - 1 e (eq) >] . 

(55) 

where S is the diagonal relaxation matrix: 

S=diag(0,-5 2 ,-5 3 ,0,-5 5 ,0,“5 7 ,-5 8 ,-5 9 ). (56) 

In simulations using the LBE method, the initial condi- 
tions provided are usually specified by velocity and pressure 
(density) fields. Often the initial condition of/ a is set to its 
equilibrium value corresponding to the given flow fields, 
with a constant density if the initial pressure field is not 
specified. The initial conditions of f a can include the first- 
order effect . The first-order effect in moment space is 
obtained through Eq. (55): 

| e < 1 >) = S" 1 MD|/ eq) ), (57) 

where D is a diagonal differential operator: 

D o/ r= V. (58) 

Equation (57) is similar to Chapman-Enskog analysis of/ < a 1) . 

For the shear flow, only the initial velocity profile is 
given. The density mode is set to be uniform initially. The 
remaining modes are initialized as the following: 


P-1. 

(59a) 

:=-2 + 3 (u 2 x + u 2 y ). 

(59b) 

e= 1 — 3(«j + «y), 

(59c) 

Qx ~~ ~ u x * 

(59d) 

Qy ~ u y * 

(59e) 
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Pxx={ul-u))-Jf(dx'*x-dy u y)’ ( 59f ) 

1 

Pxy z ='* x U y -—{d y U x + d x U y ). (59g) 

The terms in p xx and p xy involving derivatives of the veloc- 
ity field take into account viscous effects in the initial con- 
ditions. These terms are obtained through Eq. (57). The first- 
order terms in turn induce second-order contributions (with 
respect to space derivatives) which are not included here. 
This leads to weak transients of short duration if there is 
separation of time scales (2 — s%)<{2 -s 5 ). 

Our first test is the decay of a sinusoidal wave in a peri- 
odic system for various values of k. The numerical and the- 
oretical results agree with each other extremely well and 
confirm the k dependence of g and v. The agreement indi- 
cates that our local analysis is indeed sufficiently accurate in 
this case. 

The next case considered is more interesting and reveal- 
ing because the initial velocity contains shocks. Consider a 
periodic domain of size N x XN y — 84X4. At time f=0, we 
take a shear wave u y (x 7 0) of rectangular shape (discontinui- 
ties in u y at x~N x !4 and x—3N x IA): 

u y (x,0) = Uq, Kx^NJ 4, 

3N x /4<x^N x , 

u y (jt,0) = - U 0 , N x /4<x^3N x /4. 

The initial condition «*(*,()) is set to a constant everywhere. 
We consider two separate cases with and without a constant 
streaming velocity V. 

B. Steady case (V=0) 

For the case of zero streaming velocity, the initial condi- 
tion for u x is zero in the system. The solution of the Navier- 
Stokes equation for this simple problem is 

w y (x,/) = 2 a n ex P(“ ^^Ocos(£„a*), (60) 

n 

where a n is the Fourier coefficient of the initial velocity pro- 
file w y (x,0), v n ^v{k n ), and k n — 2rr{2n - 1)/A^. The mag- 
nitude of the u y (x, 0), Uq — 0.0001 in the simulations. 

Figures 4(a) and 4(b) show the decay of the rectangular 
shear wave simulated by the normal LBE scheme and the 
LBE scheme with second-order central interpolation (with 
r=0.5, where r is the ratio between advection length S x and 
grid size A*), respectively. (The detailed analysis of LBE 
schemes with various interpolations is provided in Appendix 
B.) The lines are theoretical results of Eq. (60) with v(k n ) 
obtained numerically. The times at which the profile of 
u y {x,t) (normalized by U 0 ) shown in Fig. 4 are t— 100, 200, 
...» 500. The numerical and theoretical results agree closely 
with each other. The close agreement shows the accuracy of 
the theory. In Fig. 4(b), the overshoots at early times due to 
the discontinuous initial condition are well captured by the 
analysis. This overshoot is entirely due to the strong k de- 



ar 



x 


FIG. 4. Decay of discontinuous shear wave velocity profile 
u y (xj). The lines and symbols (X) are theoretical [Eq. (60)] and 
numerical results, respectively. Only the positive half of each ve- 
locity profile is shown. LBE model (a) with no interpolation, (b) 
with the central interpolation and r = 0.5. 

pendence of v(k) caused by the interpolation. This phenom- 
ena is not necessarily connected to the Burnett effect, as 
claimed by a previous work [46]. This artifact is also com- 
monly observed in other CFD methods involving interpola- 
tions. 

Figure 5 shows the decay of w y (x,r) at one location of 



FIG. 5. Decay of discontinuous shear wave velocity u y (x,t) at a 
location close to the discontinuity x = 3N x J4. The solid lines and 
dashed lines are theoretical and numerical results, respectively. The 
LBE scheme with no interpolation does not have an overshooting, 
whereas the LBE scheme with central interpolation and r = 0.5 has. 
The time is rescaled as r~ 2 t. 
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discontinuity, x = 3N x I4=63. We tested the normal LBE 
scheme without interpolation and the LBE scheme with 
second-order central interpolation with r— 0.5, and compared 
the numerical results with theoretical ones. Again, the nu- 
merical and theoretical results agree very well with each 
other for both cases (with and without interpolation). Note 
that the time is rescaled as r~ 2 t in the figure. It should be 
pointed out that the LBE solutions of the flow differ from the 
analytic solution of the Navier-Stokes equation in both short- 
time and long-time behavior. Interpolation causes overshoot 
in the velocity at the initial stage. Even without interpolation, 
the LBE solution does not decay (exponentially) right away. 
This is due to the variation of the viscosity with k and this 
could be interpreted as the influence of the kinetic modes. (If 
we had a vanishingly small Knudsen number, then the k 
dependence would be negligible; however, all relaxation 
rates must be smaller than 2 so that higher modes can play a 
role.) This transient behavior is due to the higher-order effect 
(of velocity gradient), as discussed previously. 

C. Streaming case (V= constant) 

We also consider the case with a constant streaming in the 
initial velocity, i.e., w /C (x,0) = ^ = 0.08. This allows us to 
check the effects of the non-Galilean invariance in the sys- 
tem. With a constant streaming velocity, the solution of the 
Navier-Stokes equation is 

Uy{xj) = ^j a n QXp(-P n kl t ) C 0 S l lc n(x~Sn V x t )^ ( 61 ) 
n 

where g n ^g(k n ) is the Galilean coefficient. 

Similarly to Fig. 4, Fig. 6 shows the evolution of u y (xj) 
for the same times as in Fig. 4. The solid lines and the sym- 
bols (X) represent theoretical and numerical results, respec- 
tively. Shocks move from left to right with a constant veloc- 
ity ^=0.08. Figures 6(a), 6(b), and 6(c) show the results for 
the normal LBE scheme without interpolation, the scheme 
with second-order central interpolation, and the scheme with 
second-order upwind interpolation, respectively. In Figs. 
6(b) and 6(c), the dashed lines are the results obtained by 
setting < g n == 1 in Eq. (61). Clearly, the effect of g(k) is sig- 
nificant. For the LBE scheme with central interpolation, the 
results in Fig. 6(b) with g( k)= 1 underpredict the overshoot- 
ing at the leading edge of the shock and overpredict the 
overshooting at the trailing edge, whereas the results in Fig. 
6(c) for the LBE scheme with upwind interpolation overpre- 
dict the overshooting at the leading edge of the shock and 
underpredict the overshooting at the trailing edge. 

VIII. CONCLUSION AND DISCUSSION 

In this paper, a generalized nine-velocity LBE model 
based on the generalized LBE model of d’Humieres [27] is 
presented. The model has the maximum number of adjust- 
able parameters allowed by the discrete velocity set. The 
values of the adjustable parameters are obtained by optimiz- 
ing the hydrodynamic properties of the model through the 
linear analysis of the LBE evolution operator. The linear 
analysis also provides the generalized hydrodynamics of the 
LBE model, from which dispersion, dissipation, isotropy, 



x 




FIG. 6. Decay of discontinuous shear wave velocity profile 
u y (x,t) with a constant streaming velocity ^ = 0.08. The solid lines 
and symbols (X) are theoretical [Eq. (61)] and numerical results, 
respectively. The dashed lines in (b) and (c) are obtained by setting 
gn = 1 in Eq. (61). LBE model (a) with no interpolation, (b) with 
central interpolation and r = 0.5, (c) with upwind interpolation and 
r = 0.5. 

and stability of the model can be easily analyzed. In sum- 
mary, a systematic and general procedure by which to ana- 
lyze the LBE models is described in detail in this paper. 
Although the model studied in this paper is relatively simple, 
the proposed procedure can be readily applied to analyze 
more complicated LBE models. 

The theoretical analysis of the model is verified through 
numerical simulation of various flows. The theoretical results 
closely predict the numerical results. The stability of the 
model is also analyzed and compared with the BGK LBE 
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model. It is found that the mechanism of separate relaxations 
for the kinetic modes leads to a model which is much more 
stable than the BGK LBE model. 

The proposed model is a Galerkin type of scheme. In 
comparison with the BGK LBE model, the proposed model 
requires the transformations between the discrete velocity 
space V and the moment space M back and forth in each step 
in the evolution equation. However, the extra computational 
cost due to this transformation is only about 10-20% of the 
total computing time. Thus, the computational efficiency is 
comparable to the BGK LBE model. Our analysis also shows 
that the LBE models with interpolation schemes have enor- 
mous numerical hyperviscosities and anisotropies due to the 
interpolations. 

We also find optimal features of the proposed nine- 
velocity model: it is difficult to improve the model by simply 
adding more velocities. For instance, we found that adding 
eight more velocities (± 1,± 2) and (±2,± 1 ) would not im- 
prove the isotropy of the model. However, our analysis does 
not provide any a priori knowledge of an optimal set of 
discrete velocities. That problem can only be solved by op- 
timization of the moment problem in velocity space [24]. It 
is also worth noting that the values of all but two (a 3 and 
y 4 ) of the adjustable parameters in our model coincide with 
the corresponding parameters in the BGK LBE model. The 
main distinction between our model and the BGK LBE 
model is that our model has the freedom to allow the kinetic 
modes to relax differently, whereas in the BGK LBE model, 
all kinetic modes relax at the same rate. This mechanism 
severely affects the stability of the BGK LBE schemes, es- 
pecially when the system is strongly overrelaxed. 

It should be mentioned that the procedure we propose 
here can be applied to analyze the linear stability of spatially 
nonuniform flows, such as the Couette flow, Poiseuille flow, 
or lid-driven cavity flow. For spatially nonuniform flows, the 
lattice Boltzmann equation is linearized over a finite domain 
including boundary conditions. This leads to an eigenvalue 
problem with many more degrees of freedom than was 
needed in the analysis of this paper. Standard Arnoldi tech- 
niques [47] allow us to determine parts of the spectrum of 
the linearized collision operator, in particular to study the 
flow stability. This analysis enables us to understand the ob- 
servation that some flows are much more stable than what is 
predicted by the linear analysis of spatially uniform flows. 
For instance, in plane Couette flow with only two nodes 
along the flow direction, the only possible values of k along 
the same direction are 0 and 7 r, which are far from the value 
of k at which the bulk instability occurs. Namely, the recip- 
rocal lattice k is not large enough to accommodate the pos- 
sible unstable modes. Furthermore, in the direction perpen- 
dicular to the flow, although the reciprocal lattice k can 
accommodate unstable shear modes, the velocity gradient, 
alters the stability of the system. (It improves the stability in 
this particular case.) 

One philosophic point must be stressed. We deliberately 
did not derive the macroscopic equations corresponding to 
the LBE model in this work; instead, we only analyzed the 
generalized hydrodynamic behavior of the modes of the lin- 
earized LBE evolution operator. We argue that if the hydro- 
dynamic modes behave exactly the same way as those of the 
linearized Navier-Stokes equations, up to a certain order of 


k, provided that the Galilean invariance is also assured up to 
a certain order of k, then we can claim that the LBE model is 
indeed adequate to simulate the Navier-Stokes equations (up 
to a certain order of k). There is no distinction between the 
LBE model and the Navier-Stokes equations up to a certain 
order of k. Thus, there is no need to use the Chapman- 
Enskog analysis to obtain the macroscopic equations from 
the LBE models. On the other hand, we have also shown 
that, in the limit of k=0, these two approaches obtain the 
same results in terms of the transport coefficients and the 
Galilean coefficient. Nevertheless, it is very difficult to apply 
the Chapman-Enskog analysis to obtain the generalized hy- 
drodynamics of the LBE models, which is important to LBE 
numerical simulations of hydrodynamic systems. The stabil- 
ity result obtained by the linear analysis presented in this 
paper is very difficult for the standard Chapman-Enskog 
analysis to obtain. Therefore, the proposed procedure by 
which to analyze the LBE model indeed contains more in- 
formation and is more general than the low-order Chapman- 
Enskog analysis. Despite its generality and power, the linear 
analysis has its limitations. Because it is a local analysis, it 
does not deal with gradients. 

Our future work will extend the analysis to fully thermal 
and compressible LBE models in three-dimensional space. 
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APPENDIX A: COUPLING BETWEEN DENSITY 
AND OTHER MODES 

To consider the coupling between the density fluctuation 
<5p = p— (p) and other modes, e , e, p xx , and p xy , the equi- 
librium values of these modes are modified as to the follow- 
ing: 


c (eq) = a 2 p + y 2 ul +j]){2-p). 

(Ala) 

e (eq) = a 3 p+ y 4 (.jl+ jy)(2- p). 

(Alb) 

^f=ri0l+/;)(2-p). 

(Ale) 


(Aid) 


where (2-p) is used to linearly approximate 1/p when the 
averaged density Po~(p) = 1* With the above modifications. 
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FIG. 7. k dependence of viscosities for various models. The 
values of the adjustable parameters and the relaxation parameters 
are the same as in Fig. 1. The solid lines, dotted lines, and dashed 
lines correspond to 0=0, tt/8, and 7r/4, respectively. LBE model 
(a) with no interpolation, (b) with central interpolation, and (c) with 
upwind interpolation. 


four elements in the first column of the linearized collision 
operator C accordingly become 


Cj2 — S 2 

Cl3 = *3 


C]8“ 


\a 2 --y 2 {V^V)) , 

^3~7A(V 2 X +V 2 y ) , 


(A2a) 

(A2b) 


(A2c) 


C K =-\s^V x V y . (A2d) 


Based on the linearized collision operator with the above 
changes, the shear and the bulk viscosities at the limit of k 
— >0 are 


1 , , / 1 1 

v 0 = -(l-V-cos *)(-- 5 


(A3) 


^o = -^(2-3c, 2 )(2-5 2 ) 

1 2 s 2 

V cos cf) ^ 

- T T — — (1-3c;)(3s 2 $ 8 -252 — 4j 8 ) 

LZC s S 2 S% 

V 2 

+ ^^[^2--'8 + 2(j 2 5 8 - J 2 -5 8 )C0S- <f>] 

V 3 cos 4> 

+ — [s 2 -s s + s 2 (s s -2)cos- 4>], (A4) 

4C s S 2 S% 

The sound modes propagate with velocity V±c s (at first or- 
der in k). The Galilean coefficient up to 0(k 2 ) is 



FIG. 8. k dependence of the Galilean coefficient g for various 
models. Solid lines, dotted lines, and dashed lines correspond to 0 
= 0, 7t/ 8, and 7r/4, respectively. The middle three curves are g(k) 
for the LBE model without interpolation, the lower three for the 
LBE model with central interpolation and r=0.5, and the upper 
three for the LBE model with upwind interpolation and r = 0.5. 


k 2 

5=1 + ^[(■yg-2)(55-i 8 )(5 5 5s _3 ' y 5 -3 ' y 8 + 6 ) 

3j 5 j 8 

k 2 V 2 

4* (cos 4 6- cos 2 #)]+ — ^[(2-“5 8 )(5 8 -5 2 )sin 2 <fi 

6c;s 2 s% 

+ 2c 2 J2(s 8 — 6j 8 + 6)cos 2 (f>]. (A5) 


APPENDIX B: INTERPOLATED LBE SCHEME 

Recently, it has been proposed to use interpolation 
schemes to interpolate {f a } from a fine mesh to a coarse 
mesh in order to improve the spatial resolution calculations 
for a limited cost in total number of nodes [43,44]. Obvi- 
ously, the interpolation schemes create additional numerical 
viscosities. The Chapman-Enskog analysis shows that any 
second- or higher-order interpolation scheme does not affect 
the viscosities in the limit k— »0 on the fine mesh. A problem 
with much greater importance in practice is to calculate the 
viscosity at finite k. To our knowledge, no such analysis is 
now available in the literature. 

In the interpolated LBE schemes, the advection step is 
altered by the interpolation scheme chosen, while the colli- 
sion step remains unchanged. The advection on a fine mesh 
combined with interpolation on a coarse mesh is the recon- 
struction step on the coarse mesh. Therefore, to obtain the 
modified linearized evolution operator L, only the advection 
operation A must be changed. In what follows, we shall con- 
sider a coarse mesh with lattice constant S x , and time step 
S t . The lattice constant of a underlying fine mesh is rS x , 
with r*£ I. Effectively, the hopping velocities of particles are 
reduced by a factor of r on coarse mesh. Therefore, dimen- 
sional analysis suggests that the sound speed is reduced by a 
factor of r, and the viscosities are reduced by a factor of r 2 in 
the limit k=0. However, the dimensional analysis does not 
provide any information about the quantitative effects of in- 
terpolation when k is finite. We shall analyze the effects of 
some commonly used second-order interpolation schemes in 
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the LBE methods. For simplicity, we shall only deal with a 
uniform mesh with square grids. 

1. Central interpolation 

The reconstruction step with second-order central interpo- 
lation is given by the following formula: 

/„( r,)- '-—■/:< Tj - *„) + (! - r ! )/;( r ( ) 

(Bl) 

where/* is the post-collision value of f a , i.e.. 


2. Upwind interpolation 

The upwind direction in the LBE method is relative to the 
particle velocity e a (the characteristics) rather than the flow 
velocity u. Therefore, the interpolation stencil is static in 
time. Second-order upwind interpolation leads to 

fa( r j) ~ ( r ; — 2 <5r 0 ) + r(2 — r)/* (r ; — 8r a ) 


, (1 —r)(2 — r) ^ 

+ « fZ(rj), 


(B7) 


where Sr a is defined in Eq. (B3). Accordingly, the phase 
factors in the advection operator given by Eq. (B4) become 


/>/.+««(/). 

and 

(B2) 

(1— r)(2 — r) , r(2 — r) , r(r-l) 
A — _ ■ ■ “ H " “i 0 '* 

2 P 2 p 2 

(B8a) 

1 

a ~ * 

(B3) 

„ (1-0(2-10, x ,r(r-l)p 2 

B- 2 +r( 2 r)p + 2 

(B8b) 

The advection operator in this case becomes 


( 1 — r)(2 — r) r(2 — r) r(r-l) 

C— + + , . 

2 q 2 q- 

(B8c) 

A= diag( \ >A,C >B ,D ,AC >CB ,BD ,DA), 

where 

(B4) 

(1 — r)(2 — r) r(r-l)q 2 

D = - j - + r(2-r)q + . 

(B8d) 


r(,+ 2 ' )p +(1 - r =) + ^i». 

(B5a) 

r( ; +1) + (i-^) + r<r : l)p . 
2 p l 

(B5b) 


(B5c) 

r ( r+1 ) , ■», . r (r- 1)<? 

: 2q +(1 r ' )+ 2 ' 

(B5d) 


where p = e ik * and q — e ik y. With the new phase factors, we 
find new results at orders 1 and 2 in k. The speed of sound 
and the Galilean coefficient are multiplied by r and the vis- 
cosity coefficients are multiplied by r 2 . 

At higher order in k, dispersion effects due to lattice arise, 
leading to differences between solutions of the standard 
Navier-Stokes equations and the flows computed using the 
LBE technique. 

As in Eq. (53), we find that the advection coefficient for 
shear waves can be made isotropic to second order in k by 
choosing 


s 5 -3r 2 


(2~fs) 

(3r 2 — s 8 ) ' 


(B6) 


which improves Eq. (53), since we can choose s 8 close to 2 
while maintaining s 5 reasonably far away from 2 (between 1 
and 3/2) by taking r 2 close to 2/3. 


where p = e lk * and q~e lk y. 

Again, the third-order term (g j) in k for the shear mode is 
anisotropic unless the following relation is satisfied: 


(2~ J 8 ) 

Si r (3r 2 -3rs s + 2s 8 )' 


(B9) 


For s g and s 5 in the usual range (s 8 near 2 and s 5 between 1 
and 3/2), the preceding equation leads to a complex value of 
r. It should be pointed out that due to the commutativity of 
propagation along x and y axes, one could apply different 
interpolation formulas along each axis, according to the 
physics of flow. For instance, a large stretch of grid can be 
applied in the direction along which flow fields do not 
change much in space, whereas in the other orthogonal di- 
rection, a normal grid (without interpolation) or even a re- 
fined grid [45] can be used, so that the aspect ratio of the 
meshes is large enough to be appropriate to the flow. 

Figure 7 shows the k dependence of the normalized shear 
viscosity v{\t)l for the LBE model with and without inter- 
polation schemes. Three orientations of k are chosen: 6— 0 
(solid line), 7r/8 (dotted line), and 77/4 (dashed line). Figures 
7(a), 7(b), and 7(c) show the v(k )fv 0 for the LBE model 
with no interpolation, with second-order central interpolation 
scheme and r= 0.5, and with second-order upwind interpo- 
lation scheme and r = 0.5, respectively. It should be stressed 
that interpolation schemes do create an enormous amount of 
numerical viscosity at k — ttI1\ Both the central and the up- 
wind interpolation schemes increase the shear viscosity at k 
- 77/2 by almost two orders of magnitude, whereas without 
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interpolation, the corresponding increase for the LBE 
scheme is at most only a factor of about 2.5 (in the direction 
0=ir/8). In all cases, the viscosity displays significant an- 
isotropy at k = / nl 2. 

Similarly to Fig. 7, Fig. 8 shows the k dependence of the 
Galilean-coefficient g(k). The three curves in the middle of 
the figure corresponding to the LBE model without interpo- 
lation. The lower three curves, g(k)=£l, correspond to the 
LBE scheme with the central interpolation, and the upper 


three curves, g(k)5=I, correspond to the LBE scheme with 
the upwind interpolation. Again, interpolations have a sig- 
nificant effect on Galilean invariance. 

One common feature observed in Figs. 7 and 8 is that the 
transport coefficients of a model along the direction of 6 
= ttIS is far from those along the directions 6=0 and 9 
= tt/ 4. This is related to the fact that for the square lattice, 
the wave vector k along the direction 9— rr/8 is not a recip- 
rocal vector of the underlying lattice. 
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ABSTRACT 

In solving high Reynolds number flow problems 
with geometrical discontinuities using the method of 
lattice Boltzmann equation, the following paradox exits. 
On the one hand, the density field is required to be 
nearly constant for the nearly incompressibility 
condition to be satisfied and the pressure near the 
geometric discontinuity is linearly proportional to the 
local density field. On the other hand, the shear stress 
and the pressure are singular near the geometry 
singularity, such as sharp comers. This often results in 
undesirable, strong local spatial oscillations which 
contaminate the solution for the flow field away from 
the singular points. The recent work by Lallemand & 
Luo [7] suggests that the use of a multi-relaxation-time 
model can improve the computational stability’ and 
reduce the undesirable dispersion. However the 
difference was considered to be largely of higher order. 
In this paper we report detailed comparison and 
assessment of the performance of these two LBE 
models: single-relaxation-time (SRT) and multi- 

relaxation-time (MRT) for various flows with 
geometric and flow singularities. Computational results 
for the pressure, viscous stresses, vorticity, and flow 
velocity in regions of large gradient show that MRT 
model significantly reduced the extent of the spatial 
oscillation near the geometric singular points and 
improved the quality of the flow field at high Reynolds 
number. The difference between the solutions of the 
two models are on the leading order in such cases. 

1. BACKGROUND 

The method of lattice Boltzmann equation (LBE) 
solves the microscopic kinetic equation for particle 
distribution function/^, £ 0» where £ is the particle 


velocity, in phase space (jc, £) and time f, from which 
the macroscopic quantities (velocity and density) are 
obtained through moment integration of fix, £ t). 
Because the solution procedure is explicit, easy to 
implement, and parallelize, the LBE method has 
increasingly become an attractive alternative 
computational method for solving fluid dynamics 
problems in various systems [1-4]. The most widely 
used lattice Boltzmann model equation is the following 
single-relaxation-time LBGK model [5], 

Mx, + e a St, t + St) t) 

= --U a (x i J)-ti eq \x t ,t)) (1) 
r 

where ffix, t) and f£ q) (x, i) are the distribution 
function and the equilibrium distribution function of 
the a-th discrete velocity respectively, r is the 
dimensionless relaxation time and e a is a discrete 
velocity vector. 

The 9-velocity (or 9-bit) LBE model on the 2-D 
square lattice (Fig. 1), denoted as the D2Q9 model, has 
been widely used for simulating 2-D flows. For 
athermal fluids, the equilibrium distributions is of the 
following form [6] 

ti eq) = pv a {\+\ e a -u + 
c 

~{e a u) 2 - ~ u u\ (2) 

2 c 4 2c" 

where w a is a weighting factor and e a is a discrete 
velocity, c = Sc/St is the lattice speed, and fix and St 
are the lattice constant and the time step, respectively. 
The discrete velocities for the D2Q9 models are 
(0, 0, 0), a = 0, rest particle 

e a =< (±l,0)c,(0,±l,)c, a = 1,3, 5,7 (3) 

(±l,±l)c, a =2, 4, 6, 8. 

and the values of the weighting factor w a are 

4/9, 
w a = < 1/9, 

1/36, 
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a = 0 
a = 1,3, 5, 7 
a = 2, 4,6, 8. 
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The density and velocity can be computed as 

e-'Zf.-'Zd*’ 


P p " 

n nt • n 


{eq) 


( 5 ) 

( 6 ) 


The speed of sound of the above LBE models is c 5 = 

c! V3 and the equation of state is that of an ideal gas 
so that 

P = pc ]■ (?) 

The viscosity of the fluid is 

v={z-V i )c]a. ( 8 ) 

With this choice for the viscosity, Eq. (1) is formally a 
second order method for solving incompressible flows 
[6]. Physical and computational stability requires that 
t> 1/2. Equation (1) is usually solved in two steps: 
collision step: f* ( x it t + St) = f a (x h t) 

(9a) 

T 

streaming step: fJjXj + e a 5t, t + St) = 

f' a {X i9 t + a) % (9b) 

which is known as the LBGK scheme [4,5]. In the 
above, * denotes the post-collision values. It is noted 
that the collision step is completely local, and the 
streaming step is uniform and requires little 
computational effort, which makes Eq. (1) ideal for 
parallel implementation. The simplicity and compact 
nature of the LBGK scheme, however, necessitate the 
use of square lattices of constant spacing {Sx = Sy), 
and consequently lead to the unity of the local 

Courant-Fridrich-Levy (CFL) number, because St = 
Sx. 

In attempting to obtain solutions for high 
Reynolds number flows using the LBE method, we 
found that the solution field for (p , u, v) often exhibit 
spatial oscillations in regions of large gradient such as 
stagnation point and sharp convex comers. Especially 
near a sharp convex comer, because the pressure and 
the vorticity are singular locally, a large gradient in the 
density or pressure field exists. Since there usully is 
insufficient resolution near the comer, the large 
gradient is often accompanied by spatial oscillations. 
Depending on the geometry, such spatial oscillation 
can propagate into the flow to contaminate the 
macroscopic variables in a large region of interest. The 
spatial oscillation may adversely affect the 
computational stability and convergence rate. 

Recently, Lallemand and Luo [7] performed 
detailed analyses on the dispersion, dissipation, and 


stability characteristics of a generalized lattice 
Boltzmann equation model proposed by d’Humieres 
[8]. It was found that by the use of multiple relaxation 
times in the generalized lattice Boltzmann equations, 
better computational stability can be achieved over the 
standard LBGK scheme due to the separation of the 
relaxations of the various kinetic modes in the 
generalized lattice Boltzmann equation model— 
hereinafter referred as the multi-relaxation-time 
(MRT) model. It is also found in Ref. [7] through the 
linearized analysis on the MRT model for various 
simple flows that the MRT model gives the same 
results, to the second order accuracy, as the single- 
relaxation-time (SRT) LBGK model does. It seems 
that these two models are equivalent in the long 
wavelength limit for macroscopic variable of interest 
and the difference is a high order effect based on their 
linear analysis. Such high order differences, however, 
can be hardly detected in simple linear flows. 

Many fluid flow problems posses complicated 
geometries and mathematical singularities. Since a 
singularity often affects numerical solutions in high 
wavenumbers, it is expected that the results of the 
MRT model be noticeably different from that of the 
SRT model, at least locally near the geometric 
singularity. For convection-dominated problems, such 
local difference in the solution behavior may also lead 
to difference in the solutions over a larger scale. It is 
important to understand how the solution based on 
MRT model behaves in such flows in comparison with 
the standard LBGK model. 

The present paper reports detailed comparison and 
assessment of the performance of these two LBE 
models for various flows with geometric and flow 
singularities. A brief background on the MRT model 
will be described first. Computational results for the 
pressure, viscous stresses, vorticity, and flow velocity 
in regions of large gradient will be compared between 
the MRT and SRT models under otherwise identical 
computational and physical parameters for: 1) Stokes 
first problem; 2) steady uniform flow over a cascade of 
zero-thickness, finite length flat plates; 3) steady 
uniform flow over a cascade of finite-thickness, semi- 
infinite length plates; 4) and steady lid-driven cavity 
flow. The flow in the Stokes first problem is singular 
at t=0. The flow in a lid-driven cavity has two singular 
comers at the intersection of the moving wall and 
stationary walls in which the viscous stresses have 
non-integrable singularities. The flows over a plate and 
a step have singularities in the pressure and stresses, 
but they are weaker than in the lid-driven cavity flow. 
These flows with varying degree of singularities allow 
for a detailed assessment of the two LBE models. The 
computational results clearly demonstrate that the 
MRT model has much better behavior in flows 
involving large gradients than the SRT LBGK model. 
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2. MULTI-RELAXATION MODEL 

In Ref. [7], a new set of variables 

R = (p,e,e,j x ,q x J y ,q y ,p xx ,p xy ) T are introduced 
and R is related to the set of 
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= MF (10) 

where M_ is the 9x9 matrix transforming F to R. In the 

vector R, p is the fluid material density, e is the 
energy, sis related to the square of the energy, y* and j y 
are the momentum density (or mass flux), q x and q y 
correspond to the energy flux, and p ^ and p xy 
correspond to the diagonal and off-diagonal 
component of the viscous stress tensor. One of the 
inherent disadvantage of the standard LBGK model is 
that everything is relaxed in the same manner as given 
by Eq. (9a). In reality, because the mass and 
momentum are conserved during the molecular 
collision, there should be no relaxation for the 
conserved quantities such as mass and momentum. 
The MRT model can take this difference into 
consideration in the design of the model. In lieu of Eq. 
(9a), the collision procedure is carried out as follows, 

e = e-s 2 (e~e eq ) , (11a) 

S* = €-S z {£-€ eq ), (lib) 

( llc ) 

q y = q y -s 7 (q y -q e y 9 ), (lid) 

p'xx=Pxx~H (Pxx~P%)’ ( lle ) 

Pxy = Pxy ~ ^ 9 iPxy ~~ Pxy ) (1 10 

where * denotes the post-collision state and the 
equilibrium values were chosen to be 


e eq =-2p + 3(u 2 +v 2 ), 

(12a) 

£ eq - p-3(u 2 +V 2 ) , 

(12b) 

qV = -« , 

(12c) 

II 

1 

< 

( 1 2d) 

p 2 = u1 - y2 . 

(12e) 

Ply = 

(120 


Before the streaming step, Eq. (9b), is carried out, one 
needs to transform the post-collision values, R , back 

to F * by using 

f'=M~'R (13) 

In writing the code, Eq. (13) can be combined with Eq. 
(1 1) to obtain a single expression 

F* =F-M~ l S(R-R') (14) 

where S is the diagonal matrix: 

S = diag( 0, s 2 > ,0, s 5 ,0, s 7 , Sg , ) . (15) 

The streaming step in the MRT model is carried out 
exactly in the same manner for each component as in 
the standard LBGK model based on Eq. (9b) 

In Ref. [7], it was shown that for the MRT model 
to give the same shear viscosity as given by Eq. (8) for 
the SRT model, one needs to set 

s % =s 9 - 1/r. (16) 

It is much more flexible to chose the rest of the 
relaxation parameters: s* s 5 , and *$ 7 * In general, 
these four parameters can be chosen to be slightly 
larger than 1. In this study, we set s 2 =s 3 =s$ -s 7 =1-2 
for simplicity. Very little difference is observed in the 
flow field if a value of 1 . 1 is used for (J 2 , s 3 , s 5 , s 7 ). It 
is worth commenting here that by setting s 2 =*3 =5 5 
=s 7 = Sz = s 9 == 1/r, the SRT model is recovered. 

3. RESULTS AND DISCUSSIONS 

To effectively demonstrate the difference between 
the solutions obtained from the SRT model and the 
MRT model, we compare various macroscopic 
variables in regions of large gradient. Four cases are 
considered. They are: i) Stokes first problem; ii) steady 
uniform flow over a cascade of zero-thickness, finite 
length flat plates at Re=1000; iii) steady uniform flow 
over a cascade of finite-thickness, semi-infinite length 
plates at Re=400 based on the inlet velocity and the 
thickness of the plate; and iv) steady lid-driven cavity 
flow at Re=1000. In all the simulations, the initial 
density is set to be po= 1. The results for the density 
(and thus pressure) are presented only in terms of the 
deviation from po or some upstream reference value. 
The results of the SRT model are obtained by running 
the same MRT code with s 2 =s 3 =s 5 =s 7 = = s 9 - 

1 It. Obviously, the relative performance between 
different LBE models depends on many factors 
including the solution characteristics, types of 
variables under investigation, and grid resolution. No 
exhaustive comparison will be made. Instead, we have 
chosen reasonable grid size in all cases to contrast the 
behavior of the two models. 
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3.1 Stoke first problem 


For an infinitely long wall to move impulsively 
with a velocity U at t=0 + , the exact solution for the 
wall shear stress is given by 

T xy, W “ J====r (17) 

^7TVt 


where p is the dynamics viscosity of the fluid. Fig. 2 
shows the relative error of the LBE solutions for the 
wall shear stress, 


I r LBE - T exact I 
I 1 xy,w * xy,w I 

exact 
* or y,w 


(18) 


The results are obtained using the relaxation time 
z^l/sg-0.55. Near t=0, there is not enough spatial 

resolution for the stokes layer of the thickness ^Jvt . 
Hence as illustrated in Fig. 2 substantial oscillations 
are present near t=0. Clearly, the error in the MRT 
model is smaller than the SRT error for t< 1 00 when 
the near-wall velocity gradient is large. Eventually, the 
effect of singularity diminishes and the two solutions 
approach each other. 


3.2 Flow over a cascade of zero-thickness, 
finite length flat plate 

The schematic of the flow is shown in the insert 
of Fig. 3a. Symmetry condition is imposed at y=± H. 
The plate is placed half-way between two gridlines so 
that the standard bounce-back condition can be used to 
update the wall condition for f a ’s. A zeroth order 
extrapolation is used at the downstream exit plane for 
fa s. A constant uniform velocity condition is imposed 
at the inlet, x/L= -2. The plate length is 40 in lattice 
unit (by taking t5c=T). The relaxation time controlling 
the shear viscosity is set to be r =0.5 12 and the 
Reynolds number based on the length is Re^ULlv 
=1000. The free stream velocity is thus (7=0.1 and 
HiL = 2 so that there 80 lattices from the plate to the 
symmetry line. 

Fig. 3a compares the density deviation, pA, as a 
function ofy at xlL =0.0125, which is half grid away 
from the leading edge, based on the MRT model and 
the SRT model under otherwise identical conditions. 
Due to the singularity in the flow at the leading edge, 
it is inevitable to have large gradients in the pressure, 
stresses, and vorticity near the leading edge at high Re. 
When there is insufficient numerical resolution, an 
unphysical spatial oscillation develops near the leading 
edge. However, MRT model is seen to be much more 
effective in suppressing the spatial oscillation for p or 
p near the leading edge. Fig. 3b compares pA as a 
function of y at x/Z=0.5125 under the same condition. 
Surprisingly, the solution based on the SRT model still 


possesses a substantial level of spatial oscillations 
even in the middle of the plate for the whole cross- 
section while the solution from the MRT model has 
become sufficiently smooth. Fig. 3c shows the viscous 
normal stress, r^, normalized by pU/L , as a function 
of y at xlL = 0.5125. Similar level of oscillations is 
observed in the SRT based solution. In this work, the 
viscous stresses are obtained using the non- 
equilibrium part of the distribution function as, 



~\e a e a Sy). (19) 

Hence no finite difference is employed for the 
evaluation of the viscous stresses. Fig. 3d compares 
the dimensionless viscous shear stress, r xyy as a 
function ofy at xiL = 0.5125. Again, the oscillations in 
the SRT based solution are noticeable outside the 
viscous boundary layer. 

To develop a further, quantitative understanding 
of the performance of the two models for flow over a 
flat plate, the streamwise variation of various 
macroscopic quantities near the plate y/L =0.0125, 
which is half-grid above the plate, are also examined. 
Fig. 4a shows the variation of the pressure coefficient 


C 


p 


P~P* 
PqU 2 1 2 ' 


( 20 ) 


at y!L =0.0125 as a function of x for solutions based on 
these two models where is the pressure at the 
centerline of the inlet. It is noted that the singularity at 
x=0 resulted in oscillation in C p for about 4-5 grid 
points after the leading edge in the MRT model. 
However, C p in the SRT based solution continues to 
oscillate across the entire plate. Fig. 4b shows 
variation of the viscous normal stress, x XXJ normalized 
by pUIL y at y/L =0.0125. The superiority of the MRT 
model over the SRT model can be clearly observed in 
regions before and after the leading edge. Fig. 4c 
compares the dimensionless wall vorticity, duldy, 
normalized by U/L y between the two models. Little 
oscillation is observed for the MRT based solution 
while the SRT based solution continues to show 
oscillatory behavior up to x/L =0.4, which cannot be 
considered as the local region of the leading edge. 


3.3 Flow over a cascade of finite thickness, 
semi-finite length plates 


The insert in Fig. 5a shows the schematic of the 
flow. In this study, H/h = 4 is used. Symmetry 
conditions are imposed at the symmetry lines at y=± 
HI 2. There are 40 grids from y=-h/2 to y=h/2 so that 
there are a total of 160 lattices between the symmetry 
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lines. The free stream velocity, U, is imposed at an 
upstream section x!h~-2 or 80 to the left of the plate. 
The downstream exit location is at xih - 8 and a zeroth 
order extrapolation is employed. The Reynolds 
number, Re=hUlv(j = 400 in this study) is based on the 
inlet velocity U and the thickness of the plate h . The 
solid walls of the plate are all located half-way 
between the grids so that the bounce-back condition is 
used to handle the no-slip condition at the solid walls. 
We used a relaxation time r =0.506 that is very close 
to 0.5 in this case. There are 40 lattices (h= 40) across 
the thickness of the plate; this corresponds to U= 0.02 
for Re= 400. 

Fig. 5a shows the streamwise variation of the x- 
component velocity at yfh= 0.5125. For x>0, this 
corresponds to half a grid above the plate. Away from 
the convex comer, the approaching flow slows down 
along y/>0.5125 due to the blockage by the plate. 
Very close to the comer, the fluid element feels the 
turning of the flow direction. According to the 
classical potential flow description, the flow 
accelerates around such a convex comer. After the 
comer (x=0), the no-slip condition causes the near- 
wall flow to slow down immediately. Thus a sharp 
drop in the velocity is observed. After a short distance 
the flow separates to form a slender bubble. While the 
MRT based solution is reasonably smooth, the SRT 
based solution exhibits significant oscillations before 
and after the comer for \x\/h<0.2. 

Fig. 5b shows the variation of the pressure 
coefficient C p at >>/Z,-0.5125. The spatial oscillation in 
the MRT based solution is rather local and is of small 
amplitude. The gradual increase of the pressure after 
the comer reflects the influence of the separation 
bubble. Further downstream, the pressure starts to 
decrease linearly as one would expect for a channel 
flow. The oscillation in the SRT based solution exists 
for |x | (h ~1 and it shows clearly that the influence of 
the comer singularity is not local in the SRT based 
solution. 

Fig. 5c compares the variation of the 
dimensionless shear stress at y/*=0.5125. While the 
MRT based solution shows no sign of oscillation in 
the upstream region of the comer, the SRT based 
solution shows a rather strong oscillation in the 
viscous shear stress all the way up to x/h —0.6. Right 
after the comer, the shear stress drops sharply due the 
geometric singularity. The insufficient resolution for 
the comer singularity resulted in oscillatory behavior 
in the shear stress downstream of the comer in the 
SRT based solution. However, the MRT based 
solution again shows very little oscillation after the 
comer. Fig. 5d compares the wall vorticity between 
the two models right after the comer. No oscillation is 
observed in the MRT model, but the SRT model 
shows visible oscillation up to x/h=0. 5. 


3.4 Lid-driven cavity flow 

The insect in Fig. 6a shows the coordinate system 
for the flow inside the cavity. The first line of the grid 
in the fluid region is at a distance ASx from the wail. 
In this study, d=0.3 is used. The boundary condition 
for zte0.5 is based on that given in Ref. [9] for curved 
geometries. The height of the cavity is ///<fr=64+2A 
With t =0.52 and /te=1000, the velocity of the moving 
wall needs to be £7=0.1032. 

The velocity field is discontinuous at the two 
comers on the moving wall. Thus the flow singularity 
is stronger than in the previous two cases where the 
velocity is continuous near the convex comer. Fig. 6a 
compares the x-component velocity as a function of y 
at x/7/=0. 00464 which is on the first grid away from 
the left vertical wall. Oscillations are observed in both 
SRT based and MRT based velocity profiles due to 
insufficient resolution for the singularity. However, 
the oscillation in the MRT solution has smaller 
amplitude and is limited to a region of 5-6 girds. The 
oscillation in the SRT solution has larger amplitude 
and propagates further into the flow field. Fig. 6b 
shows the vertical component of the velocity as a 
function of y at x/H- 0.00464. Again, the SRT 
solution has a much larger amplitude and larger region 
of the oscillation. Fig. 6c compares the velocity 
profiles of the x-component of the two solutions at the 
centerline (x!H- 0.5) in the lower half of the cavity 
with a finite difference solution based on the vorticity- 
stream-fimction formulation. It is worth noting that the 
MRT based solution is noticeably more accurate than 
the SRT based solution even in regions where one 
considers far away from the singularities. 

As a final comment, by taking advantage of many 

zero elements in A/ _1 S and recognizing various 
common factors in the expressions for the vector 
A C' S(R- R* ) , the algorithm for the collision step in 

the MRT model can be coded quite efficiently. For the 
entire computation of the collision, streaming, and the 
evaluation of macroscopic variables, the code for the 
MRT model takes only about 10% more CPU time per 
time step than an SRT code does. However, this extra 
10% work is greatly compensated by the improved 
convergence of the MRT model in suppressing 
efficiently the transient oscillation associated with the 
high-frequency pressure (acoustic) waves, the much 
improved quality of the results, and the reduced 
demand for higher resolution. 

4. CONCLUSIONS 

Based on the detailed examination of the flow 
fields in various cases, it is clear that the MRT model 
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has substantial advantages over the SRT model in 
handling the geometric singularities. While the linear 
analysis given in Ref. [7] shows that the difference 
between the SRT model and MRT model mainly lies 
in the higher order term or for the high wavenumber 
components, the present simulations in more 
complicated flows demonstrate that the difference 
between the two models can be non-local. The 
difference is substantially larger than just higher order. 
The MRT model in general provides smoother 
variations of the macroscopic quantities and has much 
smaller regions of the oscillation near a singularity. 
Since the spatial oscillation is often accompanied by 
the high frequency pressure (acoustic) waves in 
transient simulations, the MRT model also offers a 
better convergence toward steady state as well. The 
MRT model is strongly recommended. The 
quantitative comparisons will obviously depend on the 
grid resolution, and the agreement between the two 
models is expected to be improved as the grid is 
refined. However the present study has clearly 
established differences, in actual computational terms, 
between the single- and multi- relaxation time models. 
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Fig. 1 A 2-D, 9-bit (or 9-speed) lattice. 
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Fig. 2 Comparison of the relative error in the evolution of the wall 
shearstress for Stokes first problem between the SRT and MRT models. 



Fig. 3a Comparison of the density profiles near the leading edge 
(x/L=0.0125) between the SRT model and MRT model at Re=1000. 



Fig. 3b Comparison of the density profiles at x/L=0.5125 between the SRT 
model and MRT model at Re=1000. 
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Fig. 3c Comparison of the viscous normal stress profiles at x/L-0.5125 
between the SRT model and MRT model at Re=1000. 



Fig. 3d Comparison of the viscous shear stress profiles at x!L= 0.5125 
between the SRT model and MRT model at Re=1000. 



Fig. 4a Comparison of pressure coefficient as a function of x at y/L— 0.0125 
between the SRT model and MRT model. 
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Fig. 6a Comparison of the velocity profiles of x-component at x A (/— 2). 



Fig. 6b Comparison of the velocity profiles of y-component at x=A (/- 2). 



Fig. 6c Comparison of the velocity profiles of jc-component atx///=0.5 in the 
lower Region of the cavity. 
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ABSTRACT 

A finite difference-based lattice Boltzmann model, employing the 2-D, 9-speed 
square (D2Q9) lattice for the compressible Euler equations, is presented. The 
model is constructed by allowing the particles to possess both kinetic and thermal 
energies. Such a lattice structure can represent both incompressible and 
compressible flow regimes. In the numerical treatment, to attain desirable 
accuracy, the total-variation-diminishing (TVD) scheme is adopted with either the 
minmod function or a second-order corrector as the flux limiter. The model can 
treat shock/expansion waves as well as contact discontinuity. Both one- and two- 
dimensional test cases are computed, and the results are compared with the exact 
as well as other reported numerical solutions, demonstrating that there is 
consistency between macroscopic and kinetic computations for the compressible 
flow. 
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1. INTRODUCTION 


The lattice Gas Automata (LG A) model was originally presented by Frisch, 
Hasslacher and Pomeaufl] as an alternative to traditional methods for simulating 
the Navier-Stokes equation. The classical LGA employs a set of Boolean variables 
to represent the particle distribution. The development of LGA was obstructed by 
two fundamental difficulties: the dependence of convection on density and the 
dependence of pressure on velocity. Furthermore, this and similar models suffer 
from statistical noise and non-Galilean invariance. The lattice Boltzmann Methods 
(LBM) developed by Qian [2], and He and Luo [3] overcome these difficulties. In 
LBM, real numbers are used to represent the ensemble-averaged particle 
distribution function. LBM has attracted much attention as a novel method for 
simulation fluid flows, including viscous flows, multiphase fluids, magneto- 
hydrodynamics, reaction-diffusion systems, and flows through porous media [4], 

In the LBM model, the velocity space is discretized into a finite number of 
discrete values, and most models are constrained to zero Mach number flows. In 
this work, we focus on the LBM simulation of inviscid compressible flows with 
discontinuities. Recently, several authors have presented LBM models for 
compressible flows based on different approaches. For example, Hu [5] proposed a 
LBM model based on the hexagonal lattice. In this model, the particles possess one 
of the three energy levels, with a total of 13-bit in the model. Note that for the 
incompressible flow, only 7-bit are required with the same lattice structure. Similar 
to Hu's model, Yan [6] presented a 17-bit model with three energy levels and three 
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speeds on a square lattice. Sun [7,8] formulated an adaptive LBM model in which 
a particle possesses two different kinds of velocity, one is the migrating velocity, 
relating to the transport (or displacement) of a particle, and the other is the phase 
velocity, related to the momentum of a particle. Palmer [9], on the other hand, 
developed a lattice Boltzmann algorithm by modeling another scalar field 
distribution for internal energy. Suffice it to say that by creating an enlarged lattice 
structure, the thermal energy and the compressibility effect,can be treated. 
However, such an approach results in a more complicated model with a large 
number of bits, meaning that the number of equations and the computing time will 
increase. 

In this paper, based on the 2-D, 9-velocity square lattice (D2Q9) model, we 
propose an LBM model for simulating compressible flows. We do not increase the 
velocity set in each lattice, but instead revise the energy definition of particles, 
assuming that besides kinetic energy, the particles also possess thermal energy. 

The specific heat ratio y can be chosen freely, and the internal energy of gas can 
vary in a wide range, so the present model can be used to simulate varying Mach 
number flows. To improve stability and accuracy, we employ the finite difference- 
based LBM, previously reported in Mei and Shyy [10], to evolve the particle 
distribution function. In section 2, we formulate a finite difference-based LBM 
model, including derivation to ascertain the equilibrium distribution function. The 
Euler equations are recovered from the present LBM to the first order in time. 
Section 3 is about the computational procedure. Section 4 presents three numerical 
tests of one dimensional flow, and the results are compared with the exact 
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solutions and numerical results based on other schemes. A two-dimensional 
numerical test of pressure pulse propagation is also presented. 

2. CONSTRUCTION OF THE LATTICE BOLTZMANN MODEL 

As illustrated in Fig. 1, the model presented in this paper is based on the D2Q9 
lattice. We redefine the energy levels e A , Eb and £c, respectively, for the rest 
particles and two kinds of particles with different speed. The particle distribution 
function at node r and time r is represented by / a ( r,t). We define the mass, 
momentum and total energy at each node as 

= ( 1 ) 

a 

( 2 ) 

a 

~pu 2 +pE = Y i (\\e a | 2 +*.)/. (3) 

where z=l, 2 for 2 dimension, a=0,...,8; e a is the velocity vector of particles(see 
Fig.l). E is the internal energy of per unit mass. 
e a =c(cos( a-l)n/4,sin(a-l)7i/4) for ot=l, 3, 5, 7 

e a = V2 c(cos( a-l)7j/4,sin(a-l)7i/4) for a=2, 4, 6, 8, 
and eo=0. 

where c is the speed parameter. We will use the following velocity moment 
tensors[l 1] to aid the derivation of the LBM. 

Yj e ca e cg = 2c 2 S tj 

or=l,3,5,7 

% e « e * = 4c2S a 

a'=2,4,6,8 
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'^ J e ai e cg e ak e am “ C ^ ij 

a= 1,3.5, 7 

Z e * e * e «* e ®» =4c “ A (;^ ~ Sc ^ijkm 


where 5 ijkm =l, if i=j=k=m, otherwise 5 ijkm =0, Ajj km =(5 ij 5 km + § lk 5 jm + 8 im 5 jk ) 
and all the odd order moment tensors equal zero. Furthermore, the particles belong 
to one of three energy levels. e a =e A for a=0; e a =£B for a=l, 3,5,7; £ a =£c for 
ot=2,4,6,8. 

The evolution of the distribution function can be written as the BGK-type 
Lattice Boltzmann equation [2,3], 

% + e«f t = -4-(/„(r,0-/”(r,0) (6) 

dt dx i rAt 

(a=0,l,...,8)where t is the single-relaxation time, f‘ q is the local equilibrium 

distribution chosen to satisfy the macroscopic gas dynamic equation. The right- 
hand side of Eq.(6) represents the simplified collision operator representing the 
changes of / due to collisions of particles. Eq.(6), after being multiplied by 1, e a 
and e a /2+e a respectively, summed up for all a=0,...,8, in addition to considering 
Eq. (1) to Eq. (3), yields 


dp | dpu, 
dt dXj 

^^(Zw.)=o 


+ pE) + ^[(L 2 + £>«/« ] = 0 


In order to recover the gas dynamic equation, besides the conservation of mass. 


momentum and energy, 
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( 10 ) 


£/r=/> 

a 

Y,fa e« = P“i ( 11 ) 

a 

z4l ‘J+e a )f?=\pu^pE (12) 


the equilibrium distribution function must satisfy the following momentum and 
energy flux conditions [5], 

Z fa e m e €t = P"i*j + P S ij (13) 

a 

Z fa (t I e a T + £ a Kf = (~ A* 2 + pE + p)u ( (14) 

a L L 


The pressure p can be obtained from the equations of state and energy of 
prefect gas. 

p = (r -i)pE (15) 

where y is the ratio of specific heat of gas. 

We assume that f‘ q has the same functional expression as that in the 
incompressible D2Q9 LBM [2,3] 


fa =D 0 p + D 3 pu 2 a = 0 

fa = A 0 p + A i + A-ipUiU j-e^e^ + A 3 pu 2 a = 13,5,7 > 

fa = B o P + B lP U i e ci + B zP U i U j e ai e 0 j + B 3 P 11 2 a = 2 , 4 , 6,8 


(16) 


Substituting Eq.(16) into Eqs.(10)-(14), we obtain ten linear algebraic equations 
for determining the coefficients A,, Bi, Dj. Considering Eq. (4), (5), and note that 
Eq. (10)-(14) are equations of power of p and «, so all terms to each order of 
p can be collected. 

Specifically, form Eq. (10) we have 
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4 A) + 4B 0 + £\> ~ 1 

2 c 2 A 2 + 4 c 2 B 2 + 4A 3 + 4B 3 + D 3 = 0 


(17) 

(18) 


From Eq. (1 1) we have 

2c 2 A, + 4c 2 5, = 1 (19) 


From Eq. (12) we have 

c 4 A2 + 2 c 2 £ b A 2 + 2 c 2 A 3 + 4 £ b A 3 + 4 c 2 B 2 + 4c 2 £ c B 2 + 4 c 2 B 3 + 4e c B 3 + £ A D 3 = ^ (20) 


2c z A 0 +4 £ B A 0 +4 c 2 S 0 + 4£ C B 0 +£’ A £> 0 = E (21) 

From Eq. (13) and consider Eq. (15) we have 

8c 4 B 2 = 1 (22) 

2c 4 Aj - 8c 4 B, = 0 (23) 

2c 2 A 3 + 4c 4 B 2 + 4c 2 B 3 = 0 (24) 

2c 2 A 0 +4c 2 B 0 ={y-l)E (25) 


From Eq. (14) and consider Eq. (15) we have 

(c 4 +2c 2 £ b )A 1 +(4c 4 +4c 2 £- c )B, = ^u 2 +yE (26) 

To determine these coefficients from Eq.(17), (21) and (25) we have 

B 0 Si £ (27) 

4(£ c -2e B +e A ) 

A,=-2S 0 +(r-l)-^ (28) 

Lc 

D 0 =1-(4A 0 +4B 0 ) (29) 

From Eqs.(19) and (26) we have 
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( 30 ) 


~u 2 +yE--c 2 -e B 
B, =— 2 


2 c"(c‘ +2f c -2f B ) 
1 


A, — — 25. + 


2c 2 


5, = 


A, = 


8c 4 

2c 4 


From Eqs.(18) and (24) we have 


*3 = 


‘‘T 


8c 2 (f c -2e B +e A ) 

1 


A - ~25 3 - 2 

4c 

D 3 =--- t -4A 3 -4 B 3 


( 31 ) 

(32) 

(33) 

(34) 

(35) 

(36) 


Here we note that if the energy levels e A , e B , Ec and c are chosen to be constants, 

then A 2 , B 2 , A 3 , B 3 . D 3 are constants also. However, A 0 , B 0 , D 0 . A,, B] are functions 
of macroscopic variables E and u. 

Choosing time step A t as the small perturbation parameter 8, we use the 
multiscale technique and Chapman-Enskog expansion 

a _ a a 2 a 

dt~dt 0 +e d^ +£ a^ + - (37) 

fa +*?+€*/*+... (38) 

The macroscopic dynamic equations of mass, momentum and energy can be 
derived from Eqs. (6)-(15) 
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(39) 


dp t dpu i 
dt dx { 


= R t +0{e 2 ) 


3/0», djpUjUj+pSy) 
dt dxj 


= R 2 +0(£ 2 ) 


^-(-^rPu 2 + pE) + ^-(^-pu 2 u i + pEu t + pu ,) = /? 3 +0(£ 2 ) 
at 2 dx i 2 


(40) 

(41) 


where 


*, 

r 2 


= 0 


(42) 

(43) 


R,=eT ^ [ ^ ( \ PU2u > + P " J ) + ^-(S/.‘"4 I*" |J +£ «> V «)1 < 44 > 

In summary, with the present lattice Boltzmann formulation, the macroscopic 
Euler equations are recovered to the first-order in time. 


3. COMPUTATIONAL PROCEDURES 

In the BGK model, the evolution equation of particle distribution function is often 
cast in the form of 

f a (x + c a At,t + At) - f a (x, t) = --(f a (x, t ) - f ' q (*,/)) (45) 

T 

Eq.(45) recovers the NS equations in the nearly incompressible flow limit. The 
most frequently adopted BGK model, with A?=l, separates the evolution process 
into two steps: collision and streaming. The resulting scheme is based on the first- 
order upwind scheme, and contains substantial numerical viscosity. For 
incompressible, viscous flow computations, the numerical viscosity is subtracted 
from the prescribed viscosity to convert the scheme to the second-order center 
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difference one. Without physical viscosity, the smearing effect of the first-order 
upwind treatment is retained because the central difference scheme, along with the 
explicit Euler time stepping, is intrinsically unstable when the viscous effect 
diminishes. To improve the solution accuracy, modem concept based on the total 
variation diminishing (TVD) can be incorporated, as will be presented below. 

In the present LBM, there are seven parameters, c, e A , Eb, Ec t. Ax, and A t. The 
ratio of specific heats y can also be chosen according to Eq. (15), which, in turn, 
directly affects the lattice model. If we let e A =£ B =£ C =0 see Eq. (12), then our 
model becomes the incompressible LBM model. It is noted that c must be greater 
than the maximum speed of the fluid flow; otherwise it can cause a negative 
distribution function. 

In the following, we first present the finite difference schemes of Eq.(6), based 
on alternative methods to treat the flux terms. There are also different strategies for 
handling the collision term, as will be presented after the convection scheme. 


(i) LBM1: The minmod method 

Utilizing the well known minmod scheme[12]the discretized form of convection 
term of Eq.(6) can be written as 


e 


a 


%,=± F 

dx t d Xj a 




where F a = ej a ,F a (I+^,J) and F a (/, J + -) is the flux of 

2 2 
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fa at the boundaries of each cell. For the advection velocity eox^, we interpolate 
the flux F a to the boundaries of each cell, using the minmod limiter to estimate the 
fluxes [12]: 

F a (I + |,/) = F c (7,7) + |minmod[A to ^F a (7,7),A^F a (7,7fl (47) 

a=l,2, 8, 

where A jwdx F a (I,J) = F a (I + l,J)~ F a (7, 7), 

A backxF a V’J) = F a ( 7,7)- F a (7 - 1, 7) , and for the advection velocity eox<0, we 
interpolate the flux F a as 

F aV = FJIJ)-^mmmod[A backx F a (I,J),A fivdx F a (I,J)] (48) 

oc=4, 5, 6 

Using the same interpolation along the y direction we have 
F a (7,7+i) = F„(7,7) + |minmod[A iflr ^F cr (7,7),A >((y F a (7,7)] (49) 

0=2, 3, 4 

Fa (FJ - = F a (7 ,7 ) - i min mod[ A backy F a (7, 7 ), A ^ F a (7,7)] (50) 

a=6,7,8, 

where A >rf> F a (7, 7) = F„ (7, 7 + 1) - F ff (7,7), 

A backy F a (l,J) = F a (7,7)-F a (7,7 -1) . 
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(Ml LB M2: Horten 's second-order flux-correction scheme 


Instead of the minmod scheme, we have also implemented a second-order flux- 
correction scheme proposed by Harten [12]. Specifically, Harten suggested the 
following flux correction: 


F « j) = f a (i, j) + i (1 + 100 „ O', j)) | C„ | (1 ~~ | I) X 

2 Ax 

x min mod[/ a (i + 1, j) - f a (/, j), f a (i, j ) - f a (/ - 1, j)] 
where co>0 is a user-adjustable parameter and 


,j) = 


\fg(i + lj)-2f a (i,j) + f a (i- 1, 7) I 


(51) 


(52) 


fa O' + !. 7) - /a O', j) | + | fa 0‘> 7> “ fa O' “ 1. 7) I 
is a shock switch. We interpolate the flux F ax to the boundaries of each cell: 


Similar to the correction of flux along x-direction, we have flux correction of y- 
direction. 


F ay O'. 7) = /„ (i, 7) + ^ (1 + Q)d (i, ;)) | c | (1 - | c |) X 

^ Ay 

x min mod[/ a (i, 7 + 1)- f a (i, j), f a (/, y) - / a (/, j - 1)] 


(54) 


where 


@ay 0'» 7 ) 


1 fa O', 7 + 1) - 2/ g Q~, 7) + / g 0‘, 7 - 1) I 
I /a O’. 7 + 1) - fa O', 7)1 + 1 fa 0, 7) ~ f a 0. j ~ 1) I 


(55) 


and interpolate the flux F ay to the boundaries of each cell: 


F a («'. 7+-) = -[F„ (i, ; + 1 ) + /% (i\ y) - 


11% 0,7 + 1) -7%, (/, 7) I 


I /„ O', 7 +!)-/« 0,7) I 


-rr(/ a (h 7 + !)-/* 0,7))] (56) 
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As to the treatment of the source term, similar to the approach of Mei and 
Shyy [10] the left-hand-side in Eq. (6) can be treated by a 3-level second-order 


scheme 


--(/r/;)r' 


--t/r 1 -(2/t" 

T 


( 57 ) 


The extrapolation for f* q ensures that the relaxation term is at the (n+l)st time 
step. Finally, the discretized from of Eq.(6) can be expressed as 




1 


T 


(58) 


4. RESULTS AND DISCUSSIONS 

In this section, three well-known one dimensional test problems are calculated. A 
system of 200x6 lattices has been used in the one-dimensional test, and a two — 
dimensional problem is also presented. 

/. l-D: Sod Test 
The initial condition is 

(Pl, ul, Pl)=( 1,0,1) -0.5<x<0 

(Pr, «*P*M0.125,0,0.1) 0<x<0.5 

The numerical results using LBM minmod method and LBM Harten scheme 
are compared with the exact solution (Fig.2 and Fig.3). The solution contains the 
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formation of a shock wave and a contact discontinuity. As expected, Harten's 
second order flux-correction scheme can resolve the discontinuity with fewer 
points while exhibiting some oscillations in the contact discontinuity region. The 
width of the sharp gradients in the LBM solutions are compared to those obtained 
by numerical methods. For the detailed information of the macroscopic solution, 
see e.g, Ref[13], To evaluate the performance quantitatively. Table 1 lists the L] 
norm error in velocity, pressure and density fields of the two tests flows, along 
with that of the macroscopic method utilizing the Roe scheme. Figure 4 illustrates 
the L) norm error versus the lattice density. As expected, the first-order accuracy is 
attained for flows with discontinuities. The present LBM solutions are very 
competitive in comparison with other well-known approaches based on the 
macroscopic formulation. 

//• l-D: Lax Test. 

The initial condition is 

0 Pu “L Rl)=(0.445,0.698,3.528) -0.5<x<0 

(Ah u R , pr)=( 0.5,0.0,0.571) 0<x<0.5 

Using LBM minmod method, Fig.5 shows the density, pressure, velocity and 
internal energy profiles of both LBM and exact solutions. Again, the results are 
reasonable, with noticeable variations in velocity and pressure profiles in the 
contact discontinuity region. Figure 6 shows the comparison with the LBM flux- 
correction solutions; sharper solution profiles are observed. However, there are 
also more pressure variations around the contact discontinuity. 
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///. 1-D: Roe Test 


The Roe test with the following initial condition: 

(Pl u L , Pl)=( 1,-1, 1.8) -0.5<x<0 

(Pr, ur, Pr)=( 1 , 1 , 1 .8) 0<x<0.5 

The numerical results and exact solutions of LBM minmod method are shown 
in Fig. 7. Similar to other macroscopic schemes, the LBM results exist error in the 
middle of the density and energy profiles. Again, as demonstrated in Fig. 8, the 
LBM flux-correction method produces solutions of similar shapes but with sharper 
gradients. As well known in the macroscopic-based method, the numerical 
solutions exhibit a bump in energy and density profiles. 

IV. 2-D: Propagation o£ a Circular Pressure Pulse 

The computational domain [0,l]x[0,l] is divided into 100x100 square lattices with 
Ax=Ay=0.01 and time step At=0.0001. The initial conditions for pressure and 
velocities are p=l, u=v=0 everywhere, with a circular pressure pulse assigned at 
the center (0.5, 0.5): 

p=l+0.5xsech(z), where z=50xfi~ j -x 50 50 ) 2 +(y. . -y 5050 ) 2 

The evolution of pressure at the center point (0.5, 0.5), with time is compared 
between the present LBM and an implicit Lagrangian method, as shown in 
Fig. 9. Good agreement is observed. The figure depicts that the pressure at the 
center drops rapidly below 1 and then recovers gradually. Fig. 10 shows the 
density, u-velocity and pressure of the solution based on the minmod method at 
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t-0.2, in half of the domain. Fig.l 1 shows the corresponding solutions based on 
Harten’s second-order flux-correction scheme. These results agree with each other, 
as well as those of Ref. 14], 


5. CONCLUDING REMARKS 

In our computations, we find that the numerical accuracy and stability are not 
sensitive to t. But the value c, Ea> £b ar *d £c roust be chosen carefully, to ensure the 
numerical stability. We also employ the second order upwind difference scheme 
for discretizing the convection term of Eq. (6), but the numerical result exists strong 
oscillation near the discontinuities. 

Here we chose e a <0 to ensure that p and E satisfy the state equation Eq.(15) to 
ensure physical realizability. Comparing with the macroscopic gas dynamic 
equation, Eq.(6) is easer, the convection term is the linear function of /«. On the 
other hand, the choice of the parameters in the compressible lattice Boltzmann 
model needs more investigation. As indicated in the model derivation and the case 
studies, there is a lack of a clear guidance to select the most appropriate 
combinations. Nevertheless, the parameters chosen here do support the numerical 
computation to maintain physical realizability. 

The theoretical derivation of the lattice Boltzmann equations and the 
prescription of the corresponding lattice structure offer a framework to develop 
numerical techniques for flow simulations. Solution accuracy and computational 
efficiency are determined by the specific schemes adopted to solve the equation. 

The present model can simulate flows over a wide range of Mach numbers and the 
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shocks are well captured. Although the model is based on square lattice for 2-D, it 
can be easily applied to hexagonal lattice and 3-D cubic lattices, as well as 
employed on curvilinear coordinates and non-uniform grids. As already pointed 
out in Section 1, in comparison to other lattice models proposed for compressible 
flow simulations, the present lattice structure, by allowing the particles to possess 
both kinetic and thermal energies, is simpler in number of discrete energy levels. 
Furthermore, in contrast to the incompressible flow model, the present numerical 
scheme benefits more from more sophisticated discretization schemes, such as 
those based on the TVD concept, which can reduce the numerical viscosity at the 
expense of added computing cost. 

As is well stated in the literature, LBM, being simple in structure, utilizing 

only linear operators, and obviously parallelizable, is an attractive approach for gas 
dynamics problems. 
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Table 1 


Sod test and Lax test numerical Li norm errors. FDLBM compared with the Roe 
schemes come from Ref. [13] 



Sod’s test f=0.1644 

Lax’s test r=0.16 


Density Velocity Pressure 

Density Velocity Pressure 

LBM1 (minmod) 

0.00859 0.02752 0.00855 

0.03051 0.01228 0.00762 

LBM2(flux-correction) 0.00617 0.02165 0.00630 

0.02241 0.00963 0.00584 

ROE 

0.00836 0.01145 0.00666 

0.02827 0.02192 0.02655 


Here LBM1 is based on the minmod method to interpolate the convection term of 
Eq. (6), and LBM2 is based on Harten's second-order flux -correction scheme. 
Total of 200 lattices in space are employed. 




